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Abstract 

An imperfectly expanded supersonic jet, invari- 
ably, radiates both broadband noise and discrete fre- 
quency sound called screech tones. Screech tones are 
known to be generated by a feedback loop driven 
by the large scale instability waves of the jet flow. 
Inside the jet plume is a quasi-periodic shock cell 
structure. The interaction of the instability waves 
and the shock cell structure, as the former propa- 
gates through the latter, is responsible for the gen- 
eration of the tones. Presently, there are formu- 
las that can predict the tone frequency fairly ac- 
curately. However, there is no known way to pre- 
dict the screech tone intensity. In this work, the 
screech phenomenon of an axisymmetric jet at low 
supersonic Mach number is reproduced by numeri- 
cal simulation. The computed mean velocity profiles 
and the shock cell pressure distribution of the jet 
are found to be in good agreement with experimen- 
tal measurements. The same is true with the sim- 
ulated screech frequency. Calculated screech tone 
intensity and directivity at selected jet Mach num- 
ber are reported in this paper. The present results 
demonstrate that numerical simulation using com- 
putational aeroacoustics methods offers not only a 
reliable way to determine the screech tone intensity 
and directivity but also an opportunity to study the 
physics and detailed mechanisms of the phenomenon 
by an entirely new approach. 

1. Introduction 

Supersonic jet noise consists of three princi- 
pal components 1 . They are the turbulent mixing 
noise, the broadband shock associated noise and the 
screech tones. Screech tones are discrete frequency 
sound. At low supersonic Mach number, the screech 
tones are associated with the axisymmetric oscilla- 
tions of the jet and radiate principally in the up- 
stream direction. It has been known since the early 
work of Powell 2 that screech tones are generated by a 
feedback loop. Recent works 1 suggest that the feed- 
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back loop is driven by the instability waves of the 
jet flow. In the plume of an imperfectly expanded 
jet is a quasi-periodic shock cell structure. Figure 
1 shows schematically the feedback loop. Near the 
nozzle lip where the jet mixing layer is thin and 
most receptive to external excitation, acoustic dis- 
turbances impinging on this area excite the instabil- 
ity waves. The excited instability waves, extracting 
energy from the mean flow, grow rapidly as they 
propagate downstream. After propagating a dis- 
tance of four to five shock cells, the instability wave 
having acquired a large enough amplitude interacts 
with the quasi-periodic shock cells in the jet plume. 
The unsteady interaction generates acoustic radia- 
tion, part of which propagates upstream outside the 
jet. Upon reaching the nozzle lip region, they ex- 
cite the mixing layer of the jet. This leads to the 
generation of new instability waves. In this way, the 
feedback loop is closed. 

At the present time, there are reliable screech tone 
frequency prediction formulas 1,3 . However, there is 
no known way to predict tone intensity and direc- 
tivity; even if it is entirely empirical. This is not 
surprising for the tone intensity is determined by 
the nonlinearities of the feedback loop. 

The principle objective of the present work is to 
simulate the screech phenomenon numerically for 
low supersonic Mach number jets. It will be shown 
that numerical simulation is an accurate method 
for screech tone intensity and directivity prediction. 
Numerical simulation of jet noise generation is not a 
straightforward undertaking. Tam 4 had earlier dis- 
cussed some of the major computational difficulties 
anticipated in such an effort. First of all, the prob- 
lem is characterized by very disparate length scales. 
For instance, the acoustic wavelength of the screech 
tone is over 20 times larger than the initial thick- 
ness of the jet mixing layer that supports the insta- 
bility waves. Further, there is also a large disparity 
between the magnitude of the velocity of the radi- 
ated sound and that of the jet flow. Typically, they 
are five to six orders different. To be able to com- 
pute accurately the instability waves and the radi- 
ated sound, a highly accurate computational aeroa- 
coustics algorithm with shock capturing capability 
as well as a set of high quality numerical boundary 
conditions are required. 
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The rest of this paper is as follows. In Section 2, 
the mathematical model, the computation algorithm 
and grid design are discussed. Section 3 describes 
the various numerical boundary conditions used in 
the simulation. Section 4 elaborates on the distribu- 
tion of artificial selective damping incorporated in 
the computation algorithm. The artificial selective 
damping terms are for the elimination of the high 
wavenumber spurious waves. They have no effect 
on the low wavenumber component (the physical so- 
lution) of the computation. They help to maintain 
a high quality numerical solution free from contam- 
ination by spurious waves and numerical instabil- 
ity. Comparisons between numerical results and ex- 
perimental measurements are provided in Section 5. 
These include the mean velocity profiles, the shock 
cell structure, the dependence of the screech tone fre- 
quency on jet Mach number and screech tone inten- 
sity. Excellent agreements with experimental mea- 
surements are found. Computed directivities for the 
first two harmonics of the dominant screech tone will 
also be provided. 

2. Mathematical Model, Computation 
Scheme and Grid Design 

In this work, we are interested in simulating the 
axisymmetric mode jet screech in the jet Mach num- 
ber range of 1.0 to 1.25. The axisymmetric mode 
is the dominant screech mode for axisymmetric jets 
from convergent nozzles at these Mach numbers. For 
this purpose, only two dimensional computation in 
the x — r plane, where (r, 0 y x) are the cylindrical 
coordinates, are necessary. 

2.1. The Mathematical Model 

For an accurate simulation of jet screech genera- 
tion, it is essential that the feedback loop be mod- 
eled and computed correctly. The important ele- 
ments that form the feedback loop are the shock cell 
structure, the large scale instability wave, and the 
feedback acoustic waves. Since turbulence in the jet 
plays only an indirect role in the feedback loop, no 
attempt is made here to resolve it computationally. 
However, turbulence in the mixing layer of the jet is 
responsible for its spreading and the spreading rate 
of the jet affects the spatial growth and decay of 
the instability wave. To ensure a good simulation 
of the spreading rate, the k — e turbulence model is 
adopted. In the computation, the modified k —e of 
Ref. [5], optimized for jet flows, are used. 

Figure 2 shows the physical domain to be simu- 
lated. We will use length scale = D (nozzle exit di- 
ameter), velocity scale = ctoo (ambient sound speed), 
time scale = — , density scale = poo (ambient gas 
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and v x are a^, 2g“- and a^D, respectively. The di- 
mensionless governing equations in Cartesian tensor 
notation are, 
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where 7 is the ratio of specific heats, v is the molec- 
ular kinematic viscosity, ko = 10 -6 and £0 = 10 -4 
are small positive numbers to prevent the division 
by zero. The model constants are taken to be 5 , 


C,, = 0.0874, 

C'l = 1.4, 

= 1.7 x IQ" 6 


<7* = 0.324, 
C e2 = 2.02, 


c e = 0.377, 
P r = 0.422, 


It is to be noted that for the range of Mach number 
and jet temperature considered the Pope and Sarkar 
corrections often added to the k — £ model 5 are not 
necessary. Outside the jet flow both k and £ are 
zero. On neglecting the molecular viscosity terms, 
the governing equations become the Euler equations. 

In this work, solutions of the above set of equa- 
tions are to be found numerically. For a given jet 
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operating condition, the solution is to provide the 
shock cell structure in the jet plume, the mean flow 
as well as the instability wave in the mixing layer 
and the acoustic field of the screech tone outside the 
jet. 

2.2. Computation scheme and grid design 

In this work, the 7-point stencil DRP scheme 4,6 is 
used to time-march the solution to a time periodic 
state corresponding to the screech cycle. The coeffi- 
cients of the scheme including those of the backward 
difference stencils are given in Ref [4] (the value of 
a 3 should be 0.0208431427703). The DRP scheme 
has proven to be nearly nondispersive over a wide 
band of wavenumbers. In the acoustic region, the 
use of 8 mesh points per wavelength would be ade- 
quate. This allows a fairly coarse grid to be used in 
the entire region outside the jet flow. 

2.3. Grid Design 

The mixing layer of the jet is very thin. The small- 
est size grid is employed here to provide needed res- 
olution. Figure 3 shows the entire computation do- 
main. The domain extends 5 diameters back from 
the nozzle exit and 35 diameters long in the in- 
direction. It is 17 diameters in the r-direction. The 
domain is divided into four blocks (or subdomains) 
as far as the grid size is concerned. In Figure 3, 
these subdomains are separated by black lines. The 
black lines represent buffer regions of 3 mesh spac- 
ings. Since acoustic waves propagate with no prefer- 
ence in direction, square grids are used. The finest 
grid in the block right downstream of the nozzle exit 
has Ax = This block is enclosed by the next 
block with Ai = ^ which, in turn, is enclosed by 
another block with Ax = y^. The outer most block 
has Ax = —. In Figure 3, the dotted curve rep- 
resents more or less the edge of the jet flow. This 
is well inside the lightly shaded region in which the 
governing equations are the k — e model turbulent 
flow equations. The full Euler equations are used in 
the unshaded region. 

The bufFer region is a narrow region around the 
boundaries of a computation subdomain of uniform 
size mesh. The change in the mesh size takes place 
in the buffer region. The basic design of the buffer 
region can be found in Ref [7] . In this work, a slightly 
improved version of the basic design is used. 

3. Numerical Boundary Conditions 

Numerical boundary conditions play a crucial role 
in the simulation of the jet screech phenomenon. Re- 
cently an in-depth review of this subject was given 


by Tam 8 . For the present problem, several types 
of numerical boundary conditions are required. In 
Figure 2, outflow boundary conditions are necessary 
along boundary AB. Along boundary BCDE y radi- 
ation condition with entrainment flow are required. 
On the nozzle wall, the solid wall boundary condi- 
tion is imposed. The jet flow is supersonic. So the 
inflow boundary condition can be prescribed at the 
nozzle exit plane. Finally, the equations of motion in 
cylindrical coordinates centered on the x-axis have 
an apparent singularity at the jet axis (r — ► 0). A 
special treatment is needed to avoid the singular- 
ity computationally. Below, a brief description of 
the different boundary conditions used in the simu- 
lations is provided. 

3.1. Radiation Boundary Conditions with 
Entrainment Flow 

For accurate numerical simulation, the numerical 
boundary conditions to be imposed along boundary 
BCDE must perform three functions. First, it must 
specify the ambient conditions for the entire compu- 
tation. This information is critical to the correct ex- 
pansion of the jet and the development of the shock 
cell structure. Second, it must allow the acoustic 
waves generated to leave the computation domain 
with minimal reflection. Third, it must generate the 
entrainment flow induced by the jet. The develop- 
ment of such a set of radiation boundary conditions 
with entrainment flow is discussed in Ref [8]. 

3.2. Outflow Boundary Conditions 

Along the outflow boundary AB, the mean flow 
is nonuniform. For this reason, the nonuniform out- 
flow boundary conditions of Tam and Dong 9 is used. 
However, as the outflow boundary is only 30 jet 
diameters downstream, the instability wave ampli- 
tude, although damped at such a far distance, re- 
mains quite large. To allow for weak nonlinearities, 
we nonlinearized the outflow boundary conditions 
of Tam and Dong by replacing the linear terms by 
their nonlinear counterparts. In cylindrical coordi- 
nates, the complete set of outflow boundary condi- 
tions used are, 


dp dp 1 (dp dp dp 
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where V{9) = u cos 6 + a( 1 — M 2 sin 2 6 ) £ , M = - . a 
is the speed of sound and ( 9 , i?) are spherical polar 
coordinates (the x-axis is the polar axis); the origin 
of R has been taken to be at the end of the potential 
core of the jet. The last two equations above are the 
nonlinearized form of the linear asymptotic k — e 
equations (without sources), p is the static pressure 
calculated by the entrainment flow model at the edge 
of the jet flow at the outflow boundary. 


3.3. Inflow Boundary Conditions 

At the nozzle exit plane, the flow variables are 
taken to be uniform corresponding to those at the 
exit of a convergent nozzle. In addition, both k and 
€ are assumed to be zero. In other words, the mixing 
layer is regarded to be very thin. Thies and Tam 5 , in 
their jet mean flow calculation work, found that this 
is a reasonably good way to initiate the computation. 
For cold jets, the mixing layer evolves rapidly into 
a quasi-similar state resembling that in a physical 
experiment. 

3.4. Boundedness Treatment at the Jet Axis 

In cylindrical coordinates, the governing equation 
has an apparent singularity at the jet axis (r — ► 0). 
Ref [8] discusses two ways to treat this problem. In 
the present work, the governing equations are not 
used at r = 0. Instead, the formal limit of these 
equations as r — + 0 are used. Our experience is that 
this can be implemented in a straightforward man- 
ner by extending the 7-point stencil into the negative 
region of r. The flow variables p t p and u in the r < 0 
region are determined by symmetric extension about 
the jet axis while v is obtained by an antisymmet- 
ric extension. These are the proper extensions for 
axisymmetric jet screech oscillations. 

3.5. Wall Boundary Conditions 

On the nozzle wall, the boundary condition of 
no through flow is implemented by the ghost point 
method of Tam and Dong 10 . For the purpose of 
eliminating the generation of spurious waves, extra 
amounts of artificial selective damping are imposed 
around the nozzle wall region. By judging from the 
computed results, this is an effective way to avoid 
the generation of short spurious waves. 

4. Artificial Selective Damping 

The DRP scheme is a central difference scheme 
and, therefore, has no intrinsic dissipation. For 


the purpose of eliminating spurious short waves and 
to improve numerical stability, artificial selective 
damping terms 11 are added to the discretized finite 
difference equations. 

In the interior region, the seven-point damping 
stencil 4 (with half-width a — 0.27t) is used. An in- 
verse mesh Reynolds number w ^ ere 

u a is the artificial kinematic viscosity) of 0.05 is pre- 
scribed over the entire computation domain. This 
is to provide general background damping to elimi- 
nate possible propagating spurious waves. Near the 
boundaries of the computation domain where a 7 
points stencil does not fit, the 5 and 3 points damp- 
ing stencils given in Ref. [4] are used. 

Spurious numerical waves are usually generated 
at the boundaries of a computation domain. The 
boundaries are also favorite sites for the occurrence 
of numerical instability. This is true for both exte- 
rior boundaries as well as internal boundaries such 
as the nozzle walls and buffer regions where there is 
a change in mesh size. To suppress both the gen- 
eration of spurious numerical waves and numerical 
instability, additional artificial selective damping is 
imposed along these boundaries. Along the inflow 
(radiation) and outflow boundaries, a distribution 
of inverse mesh Reynolds number in the form of a 
Gaussian function with a half-width of 4 mesh points 
(normal to the boundary) and a maximum value of 
0.1 right at the outermost mesh points is incorpo- 
rated into the time marching scheme. Adjacent to 
the jet axis, a similar addition of artificial selective 
damping is implemented with a maximum value of 
the inverse mesh Reynolds number at the jet axis 
set equal to 0.35. On the nozzle wall, the use of a 
maximum value of additional inverse mesh Reynolds 
number of 0.2 has been found to be very satisfactory. 

The two sharp corners of the nozzle lip and the 
transition point between the use of the outflow and 
the radiation boundary condition on the right side of 
the computation domain are locations requiring fur- 
ther additional numerical damping. This is done by 
adding a Gaussian distribution of damping around 
these special points. 

As shown in figure 3, the four computation sub- 
domains are separated by buffer regions. Here ad- 
ditional artificial selective damping is added to the 
finite difference scheme. In the supersonic region 
downstream of the nozzle exit, a shock cell structure 
develops in the jet flow. In order to provide the DRP 
scheme with shock capturing capability, the variable 
stencil Reynolds number method of Tam and Shen 12 
is adopted. The jet mixing layer in this region has 
very large velocity gradient in the radial direction. 
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Because of this, the U s tencil °f fh e variable stencil 
Reynolds number method is determined by search- 
ing over the seven-point stencil in the axial direction 
and only the two immediately adjacent mesh points 
in the radial direction. An inverse stencil Reynolds 
number distribution of the form, 



(16) 

where 



0 < x < 9 
9 < x 

G ( r ) {exp[ 0.8) 2 ], 

0 < r < 0.8 
0.8 < r 

is used in all numerical simulations. It is possible 
to show, based on the damping curve (a = 0.3tt) 


that the variable damping has minimal effect on the 
instability wave of the feedback loop. Also exten- 
sive numerical experimentations indicate that the 
method used can, indeed, capture the oscillatory 
shocks in the jet plume and that the time averaged 
shock cell structure compares favorably with exper- 
imental measurements. 

5. Numerical Results and 
Comparisons with Experiments 

We have been able, using the numerical algorithm 
described above, to reproduce the jet screech phe- 
nomenon computationally. Figure 4 shows the com- 
puted density field in the x — r-plane at one in- 
stance after the initial transient disturbances have 
propagated out of the computational domain. The 
screech feedback loop locks itself into a periodic cy- 
cle without external interference. As can be seen, 
sound waves of the screech tones are radiated out 
in a region around the fourth to fifth shock cells 
downstream of the nozzle exit. Most of the promi- 
nent features of the numerically simulated jet screech 
phenomenon are in good agreement with physical 
experiments 1 3 > 1 4 - 1 5 . 

5.1. Mean Velocity Profiles and Shock Cell 
Structure 

To demonstrate that the present numerical simula- 
tion can actually reproduce the physical experiment, 
we will first compare the mean flow velocity profile of 
the simulated jet with experimental measurements. 
For this purpose, the time averaged velocity profile 
of the axial velocity of a Mach 1.2 jet from 1 diam- 
eter downstream of the nozzle exit to 7 diameters 
downstream at one diameter interval are measured 
from the numerical simulation. They are shown as 


a function of rf = ^ in figure 5 where r 0 . 5 is 

the radial distance from the jet axis to the location 
where the axial velocity is equal to half the fully 
expanded jet velocity. Numerous experimental mea- 
surements have shown that the mean velocity profile 
when plotted as a function of r)* would nearly col- 
lapse into a single curve. The single curve is well 
represented by an error function in the form, 

— = 0.5[1 — erf(<T7/* )] (17) 

Uj 

where a is the spreading parameter. Extensive jet 
mean flow data had been measured by Lau 16 . By in- 
terpolating the data of Lau to Mach 1.2, it is found 
that experimentally a is nearly equal to 17.0. The 
empirical fit, formula (17), with a = 17.0 is also 
plotted in figure 5 (the circles). As can readily be 
seen, there is good agreement between the empiri- 
cal mean velocity profile and those of the numerical 
simulation. 

One important component of the screech feedback 
loop is the shock cell structure inside the jet plume. 
To ensure that the simulated shock cells are the same 
as those in an actual experiment, we compare the 
time averaged pressure distribution along the cen- 
terline of the simulated jet at Mach 1.2 with the ex- 
perimental measurements of Norum and Brown 17 . 
Figure 6 is a plot of the simulated and measured 
pressure distribution as a function of downstream 
distance. It is clear from this figure that the first five 
shocks of the simulation are in good agreement with 
experimental measurements both in terms of shock 
cell spacing and amplitude. Beyond the fifth shock 
cell, the agreement is less good. At this time, we are 
unable to determine the cause of the discrepancy. 
However, it is known that screech tones are gener- 
ated around the fourth shock cell. Therefore, any 
minor discrepancies downstream of the fifth shock 
cell would not invalidate our screech tone simula- 
tion. 

During a screech cycle, the shock cell is not sta- 
tionary. In the past, Westley and Woolley 13 had 
made extensive high speed stroboscopic schlieren ob- 
servations of the motion of the shock cells and the 
disturbances/instability wave in the mixing layer of 
the jet. Figure 7a is their hand sketch of the promi- 
nent features inside the jet plume. Figure 7b is the 
density field in a plane cutting through the center- 
line of the simulated jet. In comparing figures 7a and 
7b one must be aware that the lighting in schlieren 
observation gives an integrated view of the density 
field across the jet. Despite this inherent difference, 
there are remarkable similarities between the twofig- 
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ures. Not only the gross features of the large tur- 
bulence structure (in the form of toroidal vortices) 
and shock cells are alike, the detailed features of the 
curved shocks are nearly the same* Based on the 
above comparisons, it is believed that the numerical 
simulation, indeed, can reproduce all the important 
elements of the screech phenomenon. . 

5.2. Screech Tone Frequency and Intensity 

It is well known that at low supersonic jet Mach 
number, there are two axisymmetric screech modes; 
the A\ and the A 2 modes. Earlier, Norum 14 had 
compared the frequencies of the A\ and A 2 modes 
measured by a number of investigators. His com- 
parison indicates that the screech frequencies and 
the Mach number at which the transition from one 
mode to the other takes place (staging) vary slightly 
from experiment to experiment. It is generally 
agreed among experimentalists that the screech phe- 
nomenon is extremely sensitive to minor details of 
the experimental facility and jet operating condi- 
tions. 

In the present numerical simulation, both the A\ 
and A 2 axisymmetric screech modes are reproduced. 
Figure 8 shows the variation of -jj, where A is the 
acoustic wavelength of the tone, with jet Mach num- 
ber obtained by the present numerical simulation. 
Since £ = where / is the screech frequency, 

this figure essentially provides the frequency Mach 
number relationship. Plotted on this figure also are 
the measurements of Ponton and Seiner 15 . The data 
from both the numerical simulation and experiment 
fall on the same two curves, one for the A\ mode and 
the other for the A 2 mode. This suggests that the 
calculated screech frequencies are in complete agree- 
ment with experimental measurements; although the 
Mach number at which staging takes place is not the 
same. 

Ponton and Seiner 15 mounted two pressure trans- 
ducers at a radial distance of 0.642D and 0.889D, 
respectively, on the surface of the nozzle lip in their 
experiment. By means of these transducers, they 
were able to measure the intensities of the screech 
tones. Their measured values are plotted in figures 
9a and 9b. The transducer of figure 9b is closer 
to the jet axis and hence shows a higher dB level. 
Plotted on these figures also are the corresponding 
tone intensities measured in the numerical simula- 
tion. The peak levels of both physical and numerical 
experiments are nearly equal. Thus, except for the 
difference in the staging Mach number, the present 
numerical simulation is, indeed, capable of providing 
accurate screech tone intensity prediction as well. 


5.3. Directivity 

The directivity patterns of the simulated screech 
tones have been measured. Typical directivities for 
the A\ and A 2 modes are shown in figures 10 and 
11. A search over the literature fails to find directiv- 
ity measurements for the axisymmetric mode screech 
tones. Validation of these results, therefore, cannot 
be carried out at this time. 

Figure 10a shows the directivity of the A\ screech 
mode (fundamental frequency) at jet Mach number 
1.18 scaled to a distance of 65 D. The directivity of 
the second harmonic is given in figure 10b. Those 
for the A 2 mode at Mach 1.2 are shown in figures 
11a and lib. Overall, the directivity patterns of the 
A\ and A 2 modes are similar. However, there are 
differences in detailed features. For the fundamen- 
tal tone, the directivity pattern consists of two prin- 
cipal lobes. One lobe radiates upstream and form 
part of the screech feedback loop. The other radi- 
ates downstream peaked at a relatively small angle 
from the jet flow direction. This is not sound gener- 
ated by the interaction of instability wave and shock 
cells. It is Mach wave radiation generated directly 
by the instability wave as it propagates down the jet 

column 1,18 . 

The directivity pattern of the second harmonic, 
figures 10b and lib, also displays two principal 
lobes. One lobe peaks around the 90 deg. direction. 
This is the principal lobe. The sound is generated by 
the nonlinearities of the source (nonlinear instabil- 
ity wave shock cell interaction). The other lobe is in 
the upstream direction. We believe this is generated 
by the nonlinear propagation effect of the upstream 
propagating feedback acoustic waves. The screech 
tone intensity is quite high. This leads immediately 
to wave steepening and the generation of harmon- 
ics. An examination of the waveforms measured at 
upstream locations confirms that they are not sinu- 
soidal but somewhat distorted. Thus the main lobes 
of the second harmonic are of entirely different ori- 
gin. 

Concluding Remarks 

Recently, rapid advances have been made in the 
development of computational aeroacoustic meth- 
ods. In this work, we have demonstrated that it is 
now possible to perform accurate numerical simula- 
tion of the jet screech phenomenon by the use of one 
of these methods, namely, the DRP scheme with ar- 
tificial selective damping. Numerical boundary con- 
ditions are also crucial to the success of the simula- 
tions. At the present time, such numerical boundary 
conditions are available in the literature. In a pre- 
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vious review 4 , it was pointed out, unlike traditional 
computational fluid dynamics problems, numerical 
simulation of jet noise generation is subjected to the 
difficulties of large length scale disparity and the 
need to resolve the many orders of magnitude dif- 
ferences in sound and flow. This work indicates that 
these problems can be overcome by a careful design 
of the computation grid and the use of an optimized 
high order finite difference scheme. 

The present work is restricted to the low super- 
sonic Mach number range for which the screech tones 
are axisymmetric. Future plans call for the exten- 
sion of the work to three dimensions to allow the 
simulation of flapping modes at higher Mach num- 
bers. 
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Figure 1. Schematic diagram of the screech tone feedback loop. 



Figure 2. The physical domain to be simulated. 
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Figure 5. Comparison between mean velocity profiles of numeri- 
cal simulation at Mj = 1.2 . 1 x/D = 1.0, , x/D = 2.0, 

, x/D = 3.0, , x/D = 4.0, , x/D = 5.0, 

x/D = 6.0, , x/D = 7.0, and U/Uj = 0.5[l - erf(ar}*)] t 

O , o = 17 from experiment (Lau 1981). 



Figure 6. Comparison between calculated time-averaged pressure 
distribution along the centerline of a Mach 1.2 cold jet and the 

measurement of Norum and Brown (1993). simulation, O 

experiment. 


Figure 3. The computation domain in the r - x plane showing 
the different subdomain and mesh sizes. 
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nozzle shock cells 

Figure 4. Density field from numerical simulation showing the generation and 
radiation of screech tone associated with a Mach 1.13, cold supersonic jet from a 
convergent nozzle. 
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(b) 


Figure 7. Unsteady shock cell structure and large scale distur- 
bances inside the jet plume at an instant, (a) Experimental ob- 
servation by Westley and Woolley (1968). (b) Numerical simu- 
lation. 




Mi 

Figure 9. Intensity of axisymmetric screech tones at the nozzle 
exit plane, (a) r/D = 0.889, (b) r/D = 0.642 . Experiment 
(Ponton and Seiner, 1992): O A\ mode, □ Ai mode. Numerical 
simulation: • A\ mode, ■ A 2 mode. 



Figure 8. Comparison between the acoustic wavelengths of sim- 
ulated screech tones and the measurements of Ponton and Seiner 
(1992). O, □ Measurements; •, ■ simulation. 
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Figure 10. Directivity of the Ai mode screech tone at Mj — 1.18, 
r = 65£>. (a) Fundamental frequency, (b) second harmonic. 
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Abstract 

Advances in Computational Aeroacoustics (CAA) 
depend critically on the availability of accurate, 
nondispersive, least dissipative computation algo- 
rithm as well as high quality numerical boundary 
treatments. This paper focuses on the recent de- 
velopments of numerical boundary conditions. In 
a typical CAA problem, one often encounters two 
types of boundaries. Because a finite computation 
domain is used, there are external boundaries. On 
the external boundaries, boundary conditions simu- 
lating the solution outside the computation domain 
are to be imposed. Inside the computation domain, 
there may be internal boundaries. On these inter- 
nal boundaries, boundary conditions simulating the 
presence of an object or surface with specific acoustic 
characteristics are to be applied. Numerical bound- 
ary conditions, both external or internal, developed 
for simple model problems are reviewed and exam- 
ined. Numerical boundary conditions for real aeroa- 
coustic problems are also discussed through specific 
examples. The paper concludes with a description of 
some much needed research in numerical boundary 
conditions for CAA. 

1. Introduction 

A physical problem is defined mathematically 
by the governing equations and boundary condi- 
tions. When the governing equations are dis- 
cretized to be solved computationally, the result- 
ing finite difference equations are usually of higher 
order than the original partial differential equa- 
tions. This is because high order schemes are 
needed to minimize numerical dispersion, an im- 
portant requirement of Computational Aeroacous- 
tics (CAA). The use of high order schemes will 
be assumed throughout this paper. High order fi- 
nite difference equations support extraneous solu- 
tions that are not solutions of the partial differ- 
ential equations. Thus to ensure a quality solu- 
tion, a set of numerical boundary conditions must 

1 Copyright (c)l997 by C.K.W. Tam. Published by the Amer- 
ican Institute of Aeronautics and Astronautics, Inc. with 
permission. 

* Distinguished Research Professor, Department of Mathemat- 
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be specified such that not only the physical bound- 
ary conditions are faithfully reproduced but also the 
amplitude of the extraneous solutions, if generated, 
would be minimized. 

A computation domain is inevitably finite in size. 
The result is that part of the physical domain is lost 
in the numerical simulation. It is, therefore, impor- 
tant that whatever takes place in the lost domain 
should have very little influence on the solution in- 
side the computation domain. If this is not the case, 
the effects must be simulated by the boundary condi- 
tions imposed on the boundaries of the computation 
domain. For exterior aeroacoustics problems, a set 
of nonreflecting or outflow boundary conditions are 
needed at the external boundaries. The purpose of 
the nonreflecting or outflow boundary conditions is 
to allow the radiated sound waves and the convected 
vorticity and entropy waves to leave the computation 
domain smoothly without reflection. 

The main objective of this paper is to provide an 
assessment of the recent advances in the formulation 
of numerical boundary conditions for aeroacoustics 
problems. In CAA, numerical boundary conditions 
are often developed for idealized model problems. In 
practical applications, they must be modified or ex- 
tended to account for the presence of a nonuniform 
and sometimes unknown mean flow. In many cases, 
the outgoing wave amplitude is not necessarily small. 
So linear boundary conditions would need to be ad- 
justed to allow the exit of nonlinear waves. Issues of 
this kind will also be examinined and discussed in 
this paper. 

Broadly speaking, CAA boundary conditions can 
be classified into six categories. They are: 

1. Radiation boundary conditions. 

2. Outflow boundary conditions. 

3. Wall boundary conditions. 

4. Impedance boundary conditions. 

5. Radiation/outflow boundary conditions with in- 
coming acoustics or vorticity waves. 

6. Radiation boundary conditions for ducted envi- 
ronments. 

The first three categories of boundary conditions are 
also needed in standard Computational Fluid Dy- 
namics (CFD). However, owing to the presence of 
acoustic and vorticity waves, the actual boundary 
conditions used in CAA are very different from those 
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used in traditional CFD. The last three categories 
of boundary conditions appear to be unique to CAA 
problems. 

The need for the above types of boundary condi- 
tions is best illustrated by considering the two com- 
putational aeroacoustics problems shown in figures 
1 and 2. Figure 1 shows the computation domain for 
numerical simulation of jet noise generation. The jet 
flow leaves the computation domain along boundary 
AB. Here the imposition of a set of outflow bound- 
ary conditions to allow the jet flow, sound, vortic- 
ity and entropy waves to exit smoothly would be 
most appropriate. Along boundary BCDE , radia- 
tion boundary conditions are required. Along the 
nozzle wall, wall boundary conditions are necessary. 
Figure 2 shows the computation domain for numeri- 
cal simulation of fan noise radiation from a jet engine 
inlet. An important component of fan noise is gen- 
erated by the interaction of the ingested vorticity 
waves and the rotor inside the engine. To suppress 
fan noise, a standard practice is to install sound ab- 
sorbing liners on the inner surface of the engine inlet 
as shown in figure 2. These liners are represented 
mathematically by an impedance boundary condi- 
tion. Along the exterior boundary CDEF , radiation 
boundary conditions with incoming vorticity waves 
are needed for the numerical simulation. Along in- 
ternal boundary AB, radiation boundary conditions 
for ducted environment are required to simulate the 
internal propagation of acoustic duct modes inside 
the jet engine. 

The rest of this paper is as follows. In Section 
2, numerical boundary conditions developed using 
idealized flow models will be examined and com- 
pared. In Section 3, boundary conditions developed 
for more realistic aeroacoustics problems are pre- 
sented. These two sections form the main part of 
this paper. Section 4 concludes with a discussion of 
the challenges and future directions of development 
in numerical boundary conditions for CAA. 

2. Boundary Conditions Based on Idealized 

Model Problems 

Most numerical boundary conditions available in 

the literature were developed for idealized model 

problems. Idealization, in some cases, are necessary 

to make it possible for a rigorous derivation of the 

boundary conditions. From the point of view that 

boundary conditions are local relations, the use of 
local approximations to formulate first-order bound- 
ary conditions is quite justified. The development of 
numerical boundary conditions for the acoustic wave 
equation has continued for many years. A recent re- 


view was given by Givoli 1 . For numerical boundary 
conditions relevant to CAA for which the Euler or 
Navier-Stokes equations are used, brief reviews can 
be found in the articles by Tam 2 and Lele 3 * * . 

2.1 Radiation/Inflow and Outflow Boundary 
Conditions 

It is well known that in a uniform mean flow the 
linearized Euler equations support three types of dis- 
turbances. They are the acoustic waves, the vor- 
ticity waves and the entropy waves. The acoustic 
waves propagate at sound speed relative to the mean 
flow. The vorticity as well as the entropy waves are 
frozen patterns convected downstream by the mean 
flow. Because of the presence of the three types of 
wave disturbances, each having distinct propagation 
characteristics, the outgoing disturbances present at 
the inflow and outflow boundaries are very differ- 
ent. At an inflow boundary, the only outgoing dis- 
turbances are acoustic waves. At an outflow bound- 
ary, in addition to the acoustic waves, both vorticity 
and entropy waves are convected out by the mean 
flow. Due to this distinctive difference, some au- 
thors choose to separate radiation/inflow boundary 
conditions and outflow boundary conditions as two 
different types of boundary conditions. Here we will 
do so whenever clarity demands. 

There have been many proposed radiation/inflow 
and outflow boundary conditions based on totally 
different considerations. For convenience, we will 
group them into five types as follows. 

(a) Characteristics Based Boundary Condi- 
tions 

Thompson 4,5 and Poinsat & Lele 6 proposed to 
treat the problem as one-dimensional near the 
boundary of the computation domain. The coor- 
dinate in the direction normal to the boundary is 
taken as the spatial coordinate. For Euler equations 
in one dimension, a full set of characteristics can 
be easily found. Thompson, Poinsat & Lele used 
these characteristics to form boundary conditions in- 
volving only outgoing waves. However, in two- or 
three-dimensional problems, there are no true char- 
acteristics. The characteristics boundary conditions 
work well for acoustic disturbances incident nearly 
normally on the boundary. They do not give good 
results at grazing angle of incidence or when there 
is a strong mean flow tangential to the boundary. 

(b) Boundary Conditions Derived from 
Asymptotic Solutions 

Bayliss & Turkel 7 - 8 , Hagstrom k Hariharan 9 
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and Tam & Webb 10 derived radiation and outflow 
boundary conditions by means of the asymptotic 
solutions of the governing equations. In the case 
of small amplitude disturbances superimposed on a 
uniform mean flow of density po> pressure po and 
velocity «o in the x-direction, the linearized Euler 
equations in two dimensions are, 


d U dE <9F 
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The nonhomogeneous term H on the right side of (1) 
represents distributed unsteady sources. By using 
Fourier- Laplace transforms, Tam & Webb 10 showed 
that the initial value problem of (1) has asymptotic 
solutions consisting of acoustic, vorticity and en- 
tropy waves. These asymptotic solutions have the 
form 
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At boundaries where there are only outgoing 
acoustic waves, a set of radiation boundary condi- 
tions can be derived by eliminating the unknown 
function F from (2) by first taking the t (time) and 
r derivatives. The resulting radiation boundary con- 
ditions are, 
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0 + 0 (r-*). (5) 


At the outflow region, the outgoing disturbances 
consist of a combination of acoustic, vorticity and 
entropy waves, that is, a direct sum of (2), (3) and 
(4). It turns out, it is possible to eliminate the un- 
known functions F y ¥ and x> an d upon using the 
linearized momentum equations of (1), to obtain the 
following set of outflow boundary conditions. 
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where (r, 0) are the polar coordinates. V (0) — 
u 0 cos 0 + a 0 ( 1 - A/ 2 sin 2 0)i, M = a 0 = 
is the speed of sound. 

(ii) Vorticity waves 
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Extensive numerical experiments testing the accu- 
racy of (5) and (6) have been carried out. The results 
indicate that radiation boundary conditions (5) and 
outflow boundary conditions (6) are extremely ef- 
fective, provided the sources are sufficiently far from 
the boundary of the computation domain. When 
there are sources located close to the boundary, the 
quality of the numerical solution is somewhat de- 
graded. 
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(iii) Entropy waves 
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In (2) to (4), the functions F, ^ and x depend on 
the initial condition and the unsteady source distri- 
bution. 


(c) Absorbing Boundary Conditions 

A different idea to deal with exterior boundary 
conditions is to use an absorbing layer. An ab- 
sorbing layer usually consists of 10 to 20 mesh 
points in which damping terms are introduced to 
damp out the incident waves. The development 
of absorbing boundary conditions has been pur- 
sued by many investigators including Engquist &; 
Majda 11 , Higdon 12 * 13 , Kosloff & Kosloff 14 and Jiang 
& Wong 15 . 

In a more recent work, the idea of absorbing the in- 
cident wave w as extended and refined by Colonius et 
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al. into a sponge and exit zone with grid stretching 
and filtering. Their work is directly related to the 
earlier work by Rai & Moin 17 . Similar proposal but 
without grid stretching was advanced before by Is- 
raeli and Orszag 18 . A somewhat different approach 
was suggested by Ta’asan & Nark 19 . They artifi- 
cially modified the governing equations in a buffer 
zone so that the mean flow becomes supersonic in the 
outward direction. This idea was further extended 
by Hayden and Turkel 20 to the full Euler equations 
in conservation form. Most recently Freund 21 pro- 
posed a zonal approach combining the absorbing 
boundary idea and the technique of Ta’asan & Nark. 

(d) Perfectly Matched Layer 

In an absorbing layer, the addition of artificial 
damping terms to the governing equations for the 
purpose of damping out the incidence disturbances 
also can lead to substantial reflections at the inter- 
face. Berenger 22 ' 23 , in his work on computational 
electromagnetics, found that it is possible to formu- 
late an absorbing layer without reflection. Such a 
layer has come to be known as a perfectly matched 
layer (PML). It has found applications in computa- 
tional aeroacoustics, elastic wave propagation 24 and 
other areas. Hu 25 was the first to apply PML to 
acoustics problems governed by the linearized Eu- 
ler equations with uniform mean flow. He has since 
extended his work to nonuniform flow and for the 
fully nonlinear Euler equations 26 . Further applica- 
tions of PML can be found in the recent works of 
Hu and coworkers 27 ’ 28 . One great advantage of the 
PML method is that if the mean flow is uniform the 
boundary of the computation domain can be put 
very close to the acoustic sources. This sometimes 
allows the use of a small computation domain. 

Although PML has been demonstrated to perform 
exceedingly well computationally yet the PML equa- 
tions with a mean flow are unstable. Consider the 
computation of small amplitude disturbances super- 
imposed on a uniform mean flow in a computation 
domain as shown in figure 3. Let’s use Ax = Ay 
(the mesh size) as the length scale, ao (the sound 
speed) as the velocity scale, ^ as the time scale, 
Po°o (where po is the mean density) as the pres- 
sure scale. The dimensionless governing equations 
in the PML are formed by splitting the linearized 
Euler equations according to the spatial derivatives. 
An absorption term is added to each of the equations 
with spatial derivative in the direction normal to the 
layer. For example, for the PML on the right bound- 
ary of figure 3, region (1), the governing equations 
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where M x and M y are the mean flow Mach num- 
bers in the x and y directions, cr is the absorption 
coefficients. 

Suppose we look for solutions with (x,y } t) depen- 
dence in the form exp[i(ax + (3y — art)]. It is easy 
to find from (7) that the dispersion relations of the 
PML region are, 
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In the limit a -+ 0, (8) and (9) become the dispersion 
relations of the acoustic and the vorticity waves of 
the linearized Euler equations. (8) is a quadric equa- 
tion in w. It has two extra roots in addition to the 
two modified acoustic modes. For small cr, the two 
spurious roots are damped but one of the modified 
acoustic roots is unstable. For larger cr, numerical 
solutions indicate that one of the spurious roots be- 
comes unstable. In any case, the equation splitting 
procedure and the addition of an absorption term, 
both are vital to the suppression of reflections at the 
interface between the computation domain and the 
PML, inadvertently lead to instabilities. 

For small a, the roots of (8) and (9) can be found 
by perturbation. Let, 

+ CTCJj 0 ^ + <r 2 u 2^ + . . • (10) 


U/00 = + (TU)^ + <t 2 u ^ + . . . (11) 

where the roots of (8) and (9) are designated by 
a superscript ‘a 1 (for acoustic waves) and t v i (for 
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vorticity waves). Substitution of (10) and (11) into 
(8) and (9), it is straightforward to find, 

c4 a) =w+, u 0, 0 (12) 

where 
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Clearly if or has positive imaginary part, 
the mode is unstable. (13) has a simple interpreta- 
tion in the case M y = 0. In this special case, (13) 
reduces to 


terms 29 to the discretized form of (7). The design of 
the artificial selective damping stencil is such that 
there is almost no damping on the long (physical) 
waves. Thus the inclusion of these terms in the fi- 
nite difference scheme should not alter the perfectly 
matched condition for the physical waves. With arti- 
ficial damping included, the discretized form of the 
first equation of (7) according to the 7-point sten- 
cil Dispersion-Relation-Preserving (DRP) scheme 10 
is (Note: all the other equations are to be treated in 
a similar way), 
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For acoustic waves with negative phase velocity; i.e., 
au± < 0 (the group velocity can, however, be posi- 
tive) the numerator of (16) is positive, there will be 
values of 0 1 for which is purely positive imagi- 
nary. Similarly, from (15), for £ < 0 and if | > 

is also purely positive imaginary. Thus the 
PML equations in the presence of a uniform flow 
with M x ^ 0 support unstable solutions. 

In a finite difference computation the dimension- 
less wavenumbers a and j3 are restricted to the range 
of — tt to tt. Following the work of Hu 25 , we will as- 
sume a PML of width equal to 10 mesh spacings. 
For a mean flow of M x = 0.3, a value of a = 1.5 
would be quite sufficient to reduce the intensity of 
the incident acoustic waves by a factor of 10 5 . Fig- 
ure 4 shows a contour map of the growth rate of the 
most unstable wave (Im(u?) is largest) in the a — (3 
plane for such a mean flow. The maximum growth 
rate is 0.035. In carrying out numerical simulation 
over a long period of time, even a weak instability 
could be a source of trouble. It is, therefore, desir- 
able to suppress the instability. One way to suppress 
the instability and, at the same time, retain per- 
fectly matched condition at the edge of the compu- 
tation domain is to add artificial selective damping 


where R& is the artificial mesh Reynolds numbers. 
By applying Fourier transform analysis to (17) and 
(18) following Ref. [29], the damping rate intro- 
duced by the last term of (17) is 


D(a, p) = -i_ dj (e-J“ + e~ W) ( 1 9) 

Ra 3 


Figure 5 shows the contours of constant damping 
rate for R& = 1.0. The coefficient dj ' s are those 
corresponding to a — 0.37T given in the appendix of 
Ref. [2]. Figure 6 shows the combined growth and 
damping rate of figures 4 and 5 for R& = 0.46. As 
can be seen, the instability is completely suppressed. 
Note that for a PML with a width of 10 mesh spac- 
ings, waves with a wavenumber a smaller than 
cannot be excited. This band of wavenumbers lies 
within the two vertical dotted lines of figure 6. 


(e) Other Methods 

In addition to the above four types of meth- 
ods, nonreflecting boundary conditions have also 
been developed by a number of investigators us- 
ing special methodology. This includes the works 
of Giles 30 , Atkins & Casper 31 , Colonius 32 , Scott et 
a/. 33 , Kroner 34 and Roe 35 . Giles used a Fourier se- 
ries approach. His work appears to have been moti- 
vated by turbomachinery noise and flow considera- 
tion. 
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(f) Evaluation of Radiation/ Inflow and Out- 
flow Boundary Conditions 

During the last few years, there have been a 
number of papers reporting the results of evalua- 
tions of the performance and accuracy of a number 
of proposed radiation and outflow boundary condi- 
tions. Hixson et a/. 36 employed a CAA problem 
with known exact solution to evaluate the quasi- 
one-dimension al characteristic boundary conditions 
of Thompson 4 - 5 , the Fourier series boundary con- 
ditions of Giles and the asymptotic boundary con- 
ditions of Tam k Webb 10 and Bayliss k Turkel 7,8 . 
They reported that the Tam k Webb boundary con- 
ditions gave satisfactory results whereas the Thomp- 
son’s boundary conditions produced significant re- 
flections. 

Hayden k Turkel 37 reported their experience 
in using the boundary condition of a number of 
investigators 4, s, 7 , 8 ^, 10 , 33 . 34 * 35 . However, the vari- 
ous proposed boundary conditions were not imple- 
mented in the computation in an identical man- 
ner. A definitive comparison becomes impossible. 
Dong 38 , in a study of radiation boundary conditions 
for nonuniform mean flow, performed a direct com- 
parison of the results using his method and those of 
the Thompson’s and Tam k Webb’s boundary con- 
ditions. The numerical results confirm the finding of 
Hixson et a/. 36 , namely, the quasi-one-dimensional 
characteristics boundary conditions can cause sig- 
nificant reflections and inaccuracies. 

It is also worthwhile to mention that two CAA 
workshops on benchmark problems have been held 
since 1994. Some of the benchmark problems 
were designed to test radiation/inflow and outflow 
boundary conditions. In each of the workshop 
proceedings 39,40 , there is a section on comparisons 
of computed results and exact (nearly exact) solu- 
tions. They provide a measure of the quality of the 
various numerical boundary conditions used. 

2.2 Wall Boundary Conditions for High- 
Order Schemes 

In CAA, high order finite difference schemes are 
used because they have less numerical dispersion. 
However, a high order finite difference equation sup- 
port spurious solutions that have no relationship to 
the original partial differential equation. These spu- 
rious solutions are unavoidably excited at a wall. 
For aeroacoustics problems, the spurious waves are 
of two types, propagating waves with short wave 
lengths and spatially damped waves. Thus when 
an acoustic wave pulse impinges on a wall, in ad- 
dition to the reflected waves, spurious short waves 


will also be emitted in a high order finite difference 
solution. Furthermore, the spatially damped waves 
would also be generated. But they decay as they 
propagate away from the wall. Effectively they form 
a numerical boundary layer on the wall surface. 

There are two major difficulties in developing wall 
boundary conditions for high order finite difference 
schemes. First, high order finite difference equations 
require additional boundary conditions, beyond the 
physical boundary conditions of the original prob- 
lem, to define a unique solution. These additional 
boundary conditions, or the way to handle the need 
for these boundary conditions, must be found so that 
only very small amplitude spurious waves are ex- 
cited. Second, in the discretized system, each flow 
variable at either an interior or boundary mesh point 
is governed by an algebraic equation (discretized 
form of the partial differential equation). The num- 
ber of unknowns is exactly equal to the number of 
equations. Thus there will be too many equations 
and not enough unknowns if it is insisted that the 
boundary conditions at the wall mesh point are sat- 
isfied also. This is, perhaps, one of the major dif- 
ferences between partial differential equations and 
finite difference equations. 

In the literature, there is an absence of suggestions 
as how to impose wall boundary condition for high 
order schemes except for the work of Tam k Dong 41 . 
They proposed to use backward difference stencils els 
a wall is approached. This eliminates the need for 
extra boundary conditions. To provide enough un- 
knowns to enforce the physical wall boundary con- 
ditions as well as to allow the discretized govern- 
ing equations to be satisfied at mesh points on the 
wall, they suggested including ghost values at ghost 
points. Ghost points are mesh points immediately 
outside the computation domain. The number of 
ghost values to be included is equal to the num- 
ber of physical wall boundary conditions per en- 
forcement point. Tam k Dong carried out an anal- 
ysis of the problem of reflection of plane acoustic 
waves by a plane wall using the ghost point method. 
They found that the intensity of the reflected spuri- 
ous short waves is largest for normal incidence but 
is less than 0.4% of the amplitude of the incident 
wave if a resolution of 10 mesh spacings per acoustic 
wavelength is used. The thickness of the numerical 
boundary layer (defined as the distance from the wall 
at which the spurious damped numerical wave solu- 
tion drops to 0.1% of the magnitude of the incident 
wave) is a little over one mesh spacing. The ghost 
point method has since been extended by Kurbatskii 
k Tam 42 for applications to curved wall surfaces us- 
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ing Cartesian mesh. Numerical results obtained in 
a number of test cases agreed well with exact solu- 
tions. 

For acoustic wave scattering problems, Chung 
& Morris 43 proposed an Impedance Mismatched 
Method (IMM). In this method, solid bodies are re- 
placed by a new fluid medium with a large char- 
acteristic impedance, pa. When the characteristic 
impedance of the new fluid medium is infinite, it can 
be shown that the incident waves are completely re- 
flected. The advantage of the IMM method is that 
the entire computation domain including the scat- 
tering bodies can be regarded as a continuous fluid 
region making the programming exceedingly simple. 
However, unlike the ghost point method, the IMM 
cannot be used for viscous problems. 

2.3 Impedance Boundary Condition 

One of the most successful methods for suppress- 
ing fan noise radiating out the inlets of jet engines is 
to install acoustic liners inside the front part of the 
engine inlet as shown in figure 2. Mathematically, a 
liner is represented by an impedance boundary con- 
dition. The impedance, Z , is a complex quantity. If 
the time dependence is taken to be e~ tu,t then Z is 
related to the two real parameters of the liner R, the 
resistance, and X, the reactance, by 

Z = R-iX. 

Ref. [44] provides a good introduction and many 
references to the impedance of liners. In the past, 
impedance boundary condition was analyzed in the 
frequency domain. For time marching computa- 
tion, an equivalent time-domain impedance bound- 
ary condition is required. 

Presently, two entirely different approaches for de- 
veloping time-domain impedance boundary condi- 
tion are available. Both approaches have limita- 
tions. Ozyoruk & Long 45,46 , following the works 
of Sullivan 47 and Penny 48 in computational electro- 
magnetics, employed the z-transform method in im- 
plementing the impedance boundary conditions in 
the time-domain. This method provides more flexi- 
bility in fitting the frequency dependence of the re- 
sistance and reactance of the liner to experimental 
measurements. Tam & Auriault 49 used a differen- 
tial formulation of time-domain impedance bound- 
ary condition. Both methods are constrained by 
spurious numerical instability. For treatment of 
broadband noise problems, the formulation of Tam 
& Auriault is restricted by numerical instability to a 
3 parameter model. Further improvements on these 
methods are obviously desirable. 


2.4 Radiation and Outflow Boundary Condi- 
tions with Incoming Acoustic and Vorticity 
Waves 

As depicted in figure 2, there are aeroacoustics 
problems for which unsteady incoming acoustic or 
vorticity waves are an important part of the prob- 
lem. For this class of problems, the boundary condi- 
tions must allow the incoming disturbances to prop- 
agate in and the outgoing disturbances to leave the 
computation domain smoothly. There are two ways 
to treat these boundary requirements. We will refer 
to them as the nonhomogeneous boundary condi- 
tions method and the split variable method. 

The nonhomogeneous boundary conditions ap- 
proach recognizes that the computed variables are 
the direct sum of the incoming and outgoing distur- 
bances. Thus on using subscripts ‘in’ and ‘out' to 
denote the part of the flow variables associated with 
the incoming and outgoing disturbances, the outgo- 
ing disturbances can be expressed as the difference 
between the computed variables and the prescribed 
incoming disturbances; e.g., 

Pout = P — Pin • (20) 

Now at the inflow boundary, the outgoing acous- 
tic waves satisfy the radiation boundary condition 
(5). Therefore, by substitution of (20) and similar 
expressions into (5), a set of nonhomogeneous radi- 
ation boundary conditions is obtained, 



m p m 

r 1 d + *+ H 

u 

l.F(0) dr dt 2rJ 

V 


-p- 


(21) 




Pin 

r i d a , r 


w in 

1^(0) dr dt 2r 


V'm 
- Pin - 


In (21) the nonhomogeneous terms on the right 
side represent the known incoming waves. In Ref. 
[42] , the plane acoustic wave scattering problem was 
calculated numerically using (21) as the boundary 
conditions. It has been found that if the compu- 
tation is to be carried out with low spatial resolu- 
tion, then an improvement in the numerical accuracy 
is obtained if the exact finite difference solution of 
the incoming disturbances is used on the right side 
of (21). At an outflow boundary, nonhomogeneous 
outflow boundary conditions similar to (21) may be 
derived from (6). 
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Another way to generate the incoming waves is 
to divide the computation domain into an interior 
and a boundary region. In the interior region, the 
computed variables are the sum of the outgoing and 
incoming disturbances. In the boundary region (3 
mesh points for the 7-point DRP scheme), the gov- 
erning equations are either the boundary conditions 
derived from asymptotic solutions of Section 2.1(b) 
or the absorbing boundary conditions of Section 
2.1(c) or the PML equations of Section 2.1(d). The 
computed variables are the outgoing disturbances 
only. Whenever a derivative stencil extends to the 
other region, the value of the variable required can 
be obtained by using (20) and similar equations. 
Here the inflow variables are given so either p or 
Pout* whichever is appropriate can be easily found. 
In this way, the incoming disturbances are generated 
at the stencil overlapping part (overlapping with the 
boundary region) of the interior region. 

2.5 Radiation Boundary Conditions for 
Ducted Environment 

For the fan noise radiation problem illustrated in 
figure 2, when the sound waves, generated by the 
cutting of the ingested vorticity waves by the ro- 
tor, reach the opening of the jet engine inlet, part of 
them are reflected back. The reflected waves would 
be propagating in the form of duct modes if the in- 
ternal area of the engine inlet varies slowly. Unlike 
acoustic waves in the free field, duct modes are dis- 
persive. They are formed by the continuous reflec- 
tion of sound waves by the walls of the duct. Their 
propagation characteristics are very different from 
acoustic waves in free space. As a result, not all the 
radiation and outflow boundary conditions discussed 
in Section 2.1 are applicable along boundary A B of 
figure 2. 

In the Second CAA Workshop on Benchmark 
Problems 40 , several benchmark problems require the 
use of radiation boundary conditions in a ducted en- 
vironment for their solutions. For single frequency 
time periodic problems, Tam et a/. 50 developed a 
set of such radiation boundary conditions using the 
duct modes as the basis. Hu and Manthey 28 , on 
the other hand, used the PML and variable splitting 
method to form such radiation boundary conditions. 
It is necessary to point out that in a ducted environ- 
ment, the dispersion relation of the PML equations 
are not the same as those given in (8) and (9). They 
are related to the duct modes. To ensure numeri- 
cal stability, artificial selective damping is again re- 
quired in the PML. The value of the artificial mesh 
Reynolds number, R&, necessary to ensure stability 


can be found in much the same way as in Section 
2.1(d). 

3. Boundary Conditions for Real Problems 

The numerical boundary conditions discussed in 
the above section are based largely on simplified 
models. Real problems, however, are generally more 
complex. In many of these problems, numerical 
boundary conditions do not simply play a single 
role such as letting the outgoing disturbances exit 
smoothly with minimal reflections. They are to per- 
form multiple tasks. In most problems that are of 
technological significance, the mean flow is nonuni- 
form. Further, because of computer memory con- 
straint and run time limitation, the size of the com- 
putation domain is usually smaller than ideal. The 
small computer domain, forcing the boundary to be 
closer to the source or objects in the flow, puts addi- 
tional demand on the design of high quality numer- 
ical boundary conditions. There does not appear 
to have a systematic way of classifying numerical 
boundary conditions for real problems. We will il- 
lustrate, by specific examples, below how some of 
the model boundary conditions can be modified and 
extended for applications in practical CAA problems 
of current interest. 

3.1 Radiation Boundary Conditions for Sim- 
ulating Jet Noise Generation 

Let us return to the computation domain for sim- 
ulating jet noise generation in figure 1. For practical 
reasons, the size of the computation domain is typ- 
ically 30 to 40 diameters in the axial direction and 
20 to 30 diameters in the radial direction. These 
dimensions are smaller than those of the anechoic 
chambers in most physical experiments. Because of 
the proximity of the computation boundary to the 
jet flow, the boundary conditions along boundary 
BCDE are burdened with multiple tasks. Obvi- 
ously, the boundary conditions must be transparent 
to the outgoing acoustic waves radiated from the jet. 
In addition, the boundary conditions must impose 
the ambient conditions on the numerical solution. 
In other words, they specify the static conditions far 
away from the jet. Furthermore, the jet entrains a 
large volume of ambient fluid. The entrainment flow 
velocity at the computation boundary is although 
small yet not entirely negligible. For high quality nu- 
merical simulation, the boundary conditions must, 
therefore, allow the entrainment flow to enter the 
computation domain smoothly as well. 

In a recent work, Tam & Dong 51 considered the 
need to formulate a set of radiation as well as out- 
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flow boundary conditions for situations where the 
mean flow was nonuniform. They provided a gener- 
alization of the asymptotic radiation boundary con- 
ditions (5) and outflow boundary conditions (6). Let 
p, u, v and p be the weakly nonuniform mean flow 
at the boundary of the computation domain, an ap- 
propriate set of radiation boundary conditions, in 3 
dimensions, was found to be, 


m p m 

1 d u 
V(9 y r)dt v 
-P- 
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where (r, ^,x) are the cylindrical coordinates, 9 is 
the polar angle (in spherical coordinates) with the 
x-axis as the polar axis, (u, v) are the velocity com- 
ponents in the axial (x) and radial (r) directions. 
V (0, r) = u cos 0 + v sin 0 -f [a 2 — ( v cos 0 — u sin 9) 2 } 3 
and a is the speed of sound. 

In their work on numerical simulation of the gener- 
ation of axisymmetric screech tones from imperfectly 
expanded supersonic jets (see Ref. [52] for a descrip- 
tion of the jet screech phenomenon) Tam &: Shen 53 
considered a computation domain nearly identical 
to that of figure 1. They used (22) as the basis 
to develop the necessary radiation-entrainment flow 
numerical boundary conditions. It was recognized 
that the entrainment flow at the boundary of the 
computation domain would be influenced by the jet 
flow outside the computation domain. To develop an 
asymptotic entrainment flow solution Tam & Shen 
divided the jet into many evenly spaced segments 
as shown in figure 7. The jet extended beyond the 
computation domain to 60 diameters downstream. 
The mass fluxes across the boundaries of each seg- 
ment was found using empirical jet flow data. The 
difference of the mass fluxes at the two ends of each 
segment of the jet gave the amount of entrainment 
flow for the particular segment. This entrainment 
was then simulated by a point source located at the 
center of the segment. The asymptotic solution for 
a point sink located on the x-axis at x, in a com- 
pressible fluid is given by (a subscript t e i is used to 
indicate entrainment flow), 


where p 00l a<x> and 7 are the ambient gas density, 
sound speed and the ratio of specific heats. Q, the 
strength of the sink, has dimensions of pood^D 2 ; D 
is the jet diameter. On replacing (p, u,t7,p) of (22) 
by (p e , u e , v e »Pe) of (23) and by summing over the 
contributions from all the sinks, the desired radia- 
tion entrainment flow boundary conditions are ob- 
tained. 

Figure 8 shows the entrainment flow streamlines 
of a Mach 1.13 cold jet from a convergent nozzle ob- 
tained by numerical simulation. It is worthwhile to 
point out that along the right-hand boundary BC\ 
the mean flow actually flows out of the computation 
domain, exactly as observed experimentally in the 
case of a free jet. This streamline pattern would be 
very different had the entrainment flow outside the 
computation domain not been included in the sink 
flow calculation. If a cut-off were imposed at the 
right boundary of the computation domain, a recir- 
culation flow pattern would emerge. This, however, 
is inconsistent with experimental observation. 

3.2 Outflow and Jet Axis Boundary Condi- 
tions for Simulating Jet Noise Generation 

Jets are inherently unstable. The instability 
waves of jets play an important role in jet noise 
generation 52 . The instability waves, once excited at 
the nozzle lip region, grow rapidly as they propagate 
in the downstream direction. Since the jet spreads 
out in the downstream direction, it follows that the 
shear gradient and hence the instability growth rate 
decreases farther and farther downstream. Eventu- 
ally the wave would reach a location downstream 
where it becomes damped. From this point on, the 
wave amplitude decreases continuously all the way 
to the outflow boundary. In the work of Tam & 
Shen 53 the outflow boundary was located at 30 jet 
diameters downstream. At this distance, the ampli- 
tudes of the decaying instability waves (sometimes 
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referred to as large turbulence structures when there 
is less coherence) are not small. To account for the 
weak nonlinearities of the outflow disturbances, it is 
possible to nonlinearize the outflow boundary con- 
ditions (6) by replacing the linearized terms by their 
nonlinear counterpart. This yields (in cylindrical co- 
ordinates) , 


dp 

dt 
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where p is the static pressure calculated by the en- 
trainment flow model at the edge of the jet flow at 
the outflow boundary. In their jet screech tones sim- 
ulation work, Tam & Shen 53 reported that (24) pro- 
vided very satisfactory numerical results. No reflec- 
tion of any significance had been detected. 

In cylindrical coordinates, the governing equations 
have an apparent singularity at the jet axis (r — ► 0). 
For instance, the continuity equation may be written 
in the form, 

^ ^ dpu l dpw p 1 = Q (25) 

dt dr dx r d<j> r 


To handle the apparent singularity, a jet axis bound- 
ary condition may be derived by taking the formal 


limit of (25) as r 0. On noting that as r 0, 
v — y 0 while w — ► 0 faster than r, the formal limit of 
(25) is, 

dp dpv dpu 

dt dr dx 


= 0 . 


(26) 


(26), which has no apparent singularity at r = 0, is 
to be enforced at all the mesh points along the jet 
axis. 

Experience indicates that the use of (26) at r = 0 
inevitably leads to the generation of spurious short 
waves at the x-axis in a time marching simulation. 
The reason for this is simply that there is an abrupt 
change in the governing equations between the jet 
axis and the first row of mesh point off the axis. 
Such discontinuous change always leads to the radia- 
tion of short waves. For problems with axisymmetry, 
one may use the half-mesh displacement method 50 
to avoid the discontinuity. The half mesh displace- 
ment method does not involve a change in governing 


equations. It depends on the extension of the com- 
putation domain to the region r < 0 by symmetry 
and antisymmetry arguments. 

3.3 Numerical Simulation of Airframe Noise 
Generation 

During landing with the wing flaps of an aircraft 
down, the unsteady flow over the airframe is an im- 
portant source of noise. In a series of experimental 
investigation, Kendall & Ahtye 54 identified a num- 
ber of airframe noise sources; referred to as the flap 
side-edge noise, gap noise and trailing edge noise. 
One possible gap noise generation mechanism is un- 
steady flow separation around the gap between the 
wing and the flap. This possibility was investigated 
using a 2-D numerical simulation by Thies, Tam & 
Reddy 55 . For simplicity, both the wing and flap 
were approximated by flat plates as shown in fig- 
ure 9. This figure, from the numerical simulation, 
shows large unsteady separation on the suction side 
of the flap. In performing the numerical simulation, 
a relatively small computation domain was used. At 
a speed of Mach 0.15 and an angle of attack of 6 
degrees, there is a steady loading on the wing-flap 
combination. The steady loading produces a dis- 
tortion on the mean flow that extends all the way 
to the boundary of the computation domain. To 
achieve a reasonably accurate simulation, the nu- 
merical boundary conditions must not only allow the 
unsteady disturbances to leave the computation do- 
main but also account for the mean flow distortion. 
Unlike the model problems of section 2 or the work 
of Ref. [51], the difficulty here in formulating a set 
of radiation boundary conditions is that the mean 
flow is unknown a priori. 

In order to take into consideration the change in 
the mean flow at the boundary of the computation 
domain due to the presence of the wing-flap com- 
bination, one can first determine the forms of the 
asymptotic solutions of both the mean flow and the 
unsteady disturbances. This can be done by solving 
the linearized Euler equations. On using the wing 
chord L as the length scale, (incoming velocity) 
as the velocity scale, as the time scale, p^ (the 
ambient gas density) as density scale and 35 

the pressure scale, the dimensionless linearized Euler 
equations are, 


dp 

dt 


dp 
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Since only the leading term is kept in (30), (31) is 
valid to order r~ 2 for large r. On the other hand, the 
asymptotic radiation boundary conditions for acous- 
tic waves in a uniform mean flow, from (5), is 


where M is the Mach number and a is the angle of 
attack. The time independent solution of (27) can 
be found by introducing a velocity potential $(ar, y) 
defined by 
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p = M 2 p. 


Substitution of (28) into (27) gives, 
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(29) can be manipulated into the Laplace equa- 
tion by introducing a rotation and dilation of coordi- 
nates. The general solution of the Laplace equation 
can be expressed in the form of a Fourier series in po- 
lar coordinates. The lowest order nontrivial solution 
for large r is in the form of a logarithmic function. 
When rewritten in the Cartesian coordinates, it is 
found, 
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A combined asymptotic boundary conditions in 2 
dimensions that reduces to (32) for the time depen- 
dent component and (31) for the time independent 
component is, 
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On following the same reasoning, it is easy to de- 
rive a corresponding set of outflow boundary con- 
ditions suitable for use in a relatively small com- 
putation domain where weakly nonuniform two- 
dimensional mean flow is present. The equations 
are, 
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where A is an unknown constant. In the gap noise 
problem, A represents the as yet unknown loading 
on the wing-flap combination. 

It is straightforward to find by substituting (30) 
into (28), after some algebra, the following asymp- 
totic results. 
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Thies et a/. 55 implemented (33) and (34) in their 
numerical simulations of gap noise and obtained very 
satisfactory results. Figure 10 shows the sound- 
pressure- level (SPL) contours in dB from the numer- 
ical simulation. The SPL contours below the wing 
form nearly concentric circles centered at the gap 
between the wing and the flap. This indicates that 
the source of noise originates from the gap region 
in agreement with the experimental observations of 
Kendall & Ahtye 54 . 
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4. Concluding Remarks 

During the last few years, a good deal of progress 
has been made in the development of numerical 
boundary conditions for CAA. Numerical examples 
have shown that many of these boundary conditions, 
when used in conjunction with high order finite dif- 
ference schemes, are capable of providing high qual- 
ity computational results. However, a closer scrutiny 
reveals that the predominant fraction of these recent 
works is devoted primarily to radiation and outflow 
boundary conditions. Other equally important types 
of boundary conditions such as wall boundary con- 
ditions, impedance boundary conditions do not ap- 
pear to have received enough attention. The need 
for these other types of boundary conditions would 
definitely be greater in the future. For they are cru- 
cial to the application of CAA methods to fan noise, 
duct acoustics, propeller and turbomachinery noise 
problems. 

In this paper, two very important items directly 
related to numerical boundary conditions have not 
been satisfactorily discussed. The first is the dis- 
cretization and implementation of the numerical 
boundary conditions. Needless to say, the discretiza- 
tion process affects the accuracy and performance 
of a proposed boundary condition in a differential 
form. The implementation of the discretized bound- 
ary condition in relation to the time marching high 
order finite difference scheme used for the interior 
points would also have a significant impact on the 
overall accuracy and stability of the numerical solu- 
tion. The second item is error estimate. From the 
point of view of designing a computational algorithm 
for the solution of a class of aeroacoustics problems, 
a priori estimate is essential. Here order of magni- 
tude estimate is not very helpful. The real need is a 
quantitative error estimate. Most unfortunately, so 
far, very little work has been done. It is hoped that 
investigators interested in CAA would accept these 
two items as their immediate challenges. 
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Abstract 

Recently, perfectly matched layer (PML) as an ab- 
sorbing boundary condition has found widespread 
applications. The idea was first introduced by 
Berenger for electromagnetic waves computations. 
In this paper, it is shown that the PML equations 
for the linearized Euler equations support unstable 
solutions when the mean flow has a component nor- 
mal to the layer. To suppress such unstable solutions 
so as to render the PML concept useful for this class 
of problems, it is proposed that artificial selective 
damping terms be added to the discretized PML 
equations. It is demonstrated that with a proper 
choice of artificial mesh Reynolds number, the PML 
equations can be made stable. Numerical examples 
are provided to illustrate that the stabilized PML 
performs well as an absorbing boundary condition. 
In a ducted environment, the wave mode are dis- 
persive. It will be shown that the group velocity 
and phase velocity of these modes can have opposite 
signs. This results in a band of transmitted waves 
in the PML to be spatially amplifying instead of 
evanescent. Thus in a confined environment, PML 
may not be suitable as an absorbing boundary con- 
dition. 

1. Introduction 

Recently, Berenger 1,2 succeeded in formulating 
an absorbing boundary condition for computational 
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electromagnetics that has the unusual characteris- 
tic that when an outgoing disturbance impinges on 
the interface between the computation domain and 
the absorbing layer surrounding it, no wave is re- 
flected back into the computation domain. In other 
words, all the outgoing disturbances are transmit- 
ted into the absorbing layer where they are damped 
out. Such a layer has come to be known as a per- 
fectly matched layer (PML). 

Since its initial development, PML has found 
widespread applications in elastic wave propaga- 
tion 3 , computational aeroacoustics and many other 
areas. Hu 4 was the first to apply PML to aeroa- 
coustics problems governed by the linearized Euler 
equations; linearized over a uniform mean flow. He 
has since extended his work to nonuniform mean flow 
and for the fully nonlinear Euler equations 5 . Further 
applications of PML to acoustics problems including 
wavemodes in ducts can be found in the most recent 
works of Hu and coworkers 6,7 . In these references, 
examples are provided that indicate that high qual- 
ity numerical solutions could be found with PML 
used as radiation or outflow boundary conditions. 

In open unbounded domains, acoustic waves are 
nondispersive and propagate with the speed of sound 
relative to the local mean flow. Inside a duct, the 
situation is completely different. Acoustic waves are 
repeatedly reflected back by the confining walls. For 
ducts with parallel walls, the continuous reflection of 
the acoustic waves by the wall leads to the forma- 
tion of coherent wave patterns called duct modes 8,9 . 
Unlike the open domain, duct modes are disper- 
sive with phase and group velocities vary with ax- 
ial wavenumber. Because of the dispersive nature 
of the duct modes many radiation boundary con- 
ditions that work well in open domains are known 
to be inappropriate for ducted environments. For 
this reason, Tam 10 in a recent review on numerical 
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boundary conditions for computational aeroacous- 
tics, suggested that boundary condition for ducted 
environment be regarded a s a category of its own. 

There are three primary objectives in this work. 
First, we intend to show that in the presence of a 
mean flow normal to a PML, the standard PML 
equations of the linearized Euler equations support 
unstable solutions. Earlier Tam 10 had pointed out 
that the PML equations with mean flow have un- 
stable solutions. However, he did not show that the 
existence of instabilities is due to the mean flow com- 
ponent normal to the layer. The origin and char- 
acteristics of these instabilities are investigated and 
analyzed. It is interesting to mention that in his 
earliest work, Hu 4 reported that his computation 
encountered numerical instability. But by applying 
numerical filtering, he was able to obtain stable so- 
lutions. In light of our finding, we believe that what 
Hu encountered was not instability of his numerical 
scheme but that his numerical solution inadvertently 
excited the intrinsic unstable solution of the PML 
equations. Not directly related to the instability of 
the PML equations, Abarbanel and Gottlieb 11 re- 
cently analyzed the electromagnetic PML equations. 
They concluded that the equations are only weakly 
well-posed. 

Second, we will show that the instability is not 
very strong, namely, the growth rates are small. 
Also the instabilities are confined primarily to short 
waves. It is, therefore, possible to suppress the insta- 
bilities by the addition of artificial selective damping 
terms 12 to the discretized PML equations. It is im- 
portant to point out that artificial selective damping 
eliminates mainly the short waves and has negligible 
effect on the long or the physical waves. Thus the 
addition of these damping terms does not effect the 
perfectly matched conditions of the PML. 

Third, we will show that a perfectly matched layer 
may not be suitable as an absorbing boundary con- 
dition for waves in a ducted flow environment. The 
major difference between acoustic waves in an open 
domain and acoustic waves inside a duct is that in an 
unbounded region acoustic waves are nondispersive 
whereas duct modes are dispersive. It will be shown 
that in the presence of a mean flow the group and 
phase velocity of the duct modes can have opposite 
signs. Because of this, a band of transmitted waves 
will actually grow spatially instead of being damped 
in the PML. In other words, the PML equations do 
not damp these wave modes as absorbing boundary 
condition ought to do. The exception is when there 
is no mean flow in the duct. In this special case, all 
the transmitted waves are spatially damped. 


In section 2, the use of PML for open domain prob- 
lems is discussed. The stability of the PML gov- 
erning equations are investigated. It will be shown 
that the addition of damping terms to form the PML 
equations can actually cause the vorticity and acous- 
tic wave modes to become unstable. The splitting 
of the variables in formulating the PML equations 
leads to a higher order system of equations. This 
higher system supports extra solutions. These extra 
or spurious solutions are found to become unstable 
when the damping coefficient is large. Numerical ex- 
amples are provided to illustrate the spread of the 
unstable solution from the PML back into the inte- 
rior of the computation domain. 

In section 3, the effect of the addition of artifi- 
cial selective damping terms to the discretized PML 
equations is investigated. It is shown that with an 
appropriate choice of mesh Reynolds number, the 
unstable solutions of the PML equations can be sup- 
pressed. Numerical examples are given to demon- 
strate the effectiveness of the modified PML as a 
radiation/outflow boundary condition. 

Section 4 deals with the theory and application of 
PML to ducted internal flow problems. An eigen- 
value analysis is carried out to show the existence 
of a band of frequency for which the PML exerts no 
damping on the acoustic duct modes. These wave 
modes actually would grow in amplitude as they 
propagate through the PML. Numerical results are 
provided to illustrate the existence of this kind of 
amplifying ducted acoustic modes. 

2. Open Domain Problems 

Let us consider the use of PML as absorbing 
boundary condition for the solution of the linearized 
Euler equations (linearized over a uniform mean 
flow) in a two-dimensional open domain as shown 
in figure 1. We will use Ax = Ay (the mesh size) as 
the length scale, ao (the sound speed) as the veloc- 
ity scale, as the time scale, pqGq (where po is the 
mean density) as the pressure scale. The dimension- 
less governing equations in the PML are formed by 
splitting the linearized Euler equations according to 
the spatial derivatives. An absorption term is added 
to each of the equations with spatial derivative in 
the direction normal to the layer. For example, for 
the PML on the right boundary of figure 1, the gov- 
erning equations are 4 , 

du\ d , x 

— +o-«i + M r ^(ui+u 2 ) 

+ -^(PI+P2) = 0 
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du 2 

~di 


+ M y ^( Ul "t* U l) — 0 


diq d t v * 

— 4- <ri>i 4- M r — (vi + V2) = o 


into the entire jP plane except for the slit ADC. But 
since /? 2 is real and positive, for subsonic mean flow 
the point /? 2 lies outside the image. Thus no value of 
w in the upper-half w-plane would satisfy equation 
(2) indicating that there is no unstable solution. 


^ + M y -^(v 1 + v 2 ) + -^( Pl + P2 ) = 0 

(1) 

dp\ d , v 

ir + ^i + M x ^( Pl + P2 ) 

+ ^(«i + « 2 ) = 0 

^+MyA( Pl+ p 2 )+A(v 1+ v 2 )=0 

where M x and M y are the mean flow Mach num- 
bers in the x and y directions, a is the absorption 
coefficient. 

Suppose we look for solutions with (z, y,t) depen- 
dence in the form exp[i(<*£ 4- fiy — cut)]. It is easy 
to find from (1) that the dispersion relations of the 
PML equations are, 


L ocM x 
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In the limit a — + 0, (2) and (3) become the well- 
known dispersion relations of the acoustic and the 
vorticity waves of the linearized Euler equations. 

2.1. Mean Flow Parallel to PML 


2.2. Unstable Solutions of the PML Equa- 
tions 

For M x 7^ 0, the PML equations support unstable 
solutions. It is to be noted that, unlike the origi- 
nal dispersion relation of the acoustic waves, equa- 
tion (2) is a quadric equation in w. It has two extra 
roots in addition to the two modified acoustic modes. 
For small cr, the two spurious roots are damped but 
one of the modified acoustic roots is unstable. For 
larger <r, numerical solutions indicate that one of the 
spurious roots becomes unstable. In any case, the 
equation splitting procedure and the addition of an 
absorption term, both are vital to the suppression 
of reflections at the interface between the compu- 
tation domain and the PML, inadvertently, lead to 
instabilities. 

For small cr, the roots of (2) and (3) can be found 
by perturbation. Let, 

u/ a ) = w^ + a w[ a ^ 4- + . . . (5) 

w^ = w^ 4- crw^ + <t 2 u>2^ 4- . . . (6) 

where the roots of (2) and (3) are designated by 

a superscript ‘a’ (for acoustic waves) and *v 9 (for 
vorticity waves). Substitution of (5) and (6) into (2) 
and (3), it is straightforward to find, 

J 0 a) = u + , w_, 0, 0 (7) 


Dispersion relations (2) and (3) behave very differ- 
ently depending on whether there is any mean flow 
normal to the PML. When the mean flow is parallel 
to the layer; i.e., M x = 0, the solutions are stable. 
This is easy to see from (3) for the vorticity wave. 
Physically, if the mean flow is parallel to the PML, 
the vorticity waves in the computation domain, be- 
ing convected by the mean flow, cannot enter the 
layer and hence would not lead to unstable solution. 

To show that for M x = 0 all the solutions of (2) 
are stable, a simple mapping will suffice. Rewrite 
(2) in the form 


where 

w± = (aM x + 0M y ) ± (a 2 + (3 2 )± (8) 
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Figure (2) shows the image of the upper-half w-plane 
in the F plane. The upper-half w-plane is mapped 


Clearly ifwj a l or w,*^ has a positive imaginary part, 
the mode is unstable. It is easy to show, especially 
in the case M y = 0, that there are always values of 
a and /? such that w^ of (10) is purely imaginary 
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and positive. Similarly, from (11) for £ < 0 and 
if I > 1S a ^ s0 P ure *y positive imaginary. 

Thus the PML equations in the presence of a uniform 
flow with M x ^ 0 support unstable solutions. 

The unstable solutions of dispersion relations (2) 
and (3) can also be found numerically. For a given 
(a,/?) the growth rates, Wj, of the unstable solutions 
can be calculated in a straightforward manner. Fig- 
ure (3) shows the Wi contours of the most unstable 
solution of equation (2), the acoustic mode, in the 
a — /?-plane for the case M x = 0.3, M y = 0.0 and 
a — 1.5. Figure 4 shows a similar plot for the vortic- 
ity wave mode (equation (3)). In these figures only 
the unstable regions are shown. It is clear that there 
are instability waves over a wide range of wavenum- 
bers. Numerical results indicate that, in general, the 
unstable regions expand as the flow Mach number or 
the damping coefficient a increases. 


the excited unstable solution in the PML. Finally, 
figure 5d (at t = 130) shows the spread of the un- 
stable solution back into the interior computation 
domain. Figure 6 gives the corresponding waveform 
of the vorticity wave pulse. Figure 6d clearly indi- 
cates that the spread of the unstable vorticity waves 
in the PML can quickly contaminate the entire com- 
putation domain. 

Figures 7 and 8 are similar plots illustrating the 
excitation of the acoustic mode unstable solution in 
the PML. The Mach number and damping coeffi- 
cient are M x = 0.5, M y = 0.0 and cr — 1.5. The 
initial disturbance consists of a pressure pulse given 

by, 


p — p — exp 



(ii) 


u = v = 0. 


2.3. Numerical Examples 

The nature and characteristics of the unstable 
waves associated with the acoustic mode and the 
vorticity mode are quite different. To illustrate the 
excitation of these unstable solutions in the PML by 
disturbances propagating or convecting from the in- 
terior computation domain, a series of numerical ex- 
periments has been carried out. Figure 5 shows the 
results of the case of a vorticity pulse convected into 
the PML when M x — 0.3, M y — 0.2 and a — 1.0. 
The initial conditions for the pulse are (same as the 
initial conditions used by Tam & Webb 13 ) 


P = P = o 


u = 0.04yexp 



( 10 ) 


v = -0.04xexp 



The DRP time marching scheme 13 is used in the 
simulation. The PML region extends from x = 20 
to the right boundary of the computation domain. 
At the outermost boundary, the boundary condition 

p x = p 2 = p\ — P2 = Ui = «2 = v l = v 2 — 0 are 

imposed. Plotted in figure 5 are contours of the u 
velocity component. Figure 5a shows the initial pro- 
file of the contours at t = 0. Figure 5b, at t = 50, 
reveals that there is damping of the vorticity pulse 
as it begins to enter the PML. This damping is the 
result of the built-in damping, cr, of the PML. Fig- 
ure 5c, at a later time t = 90, shows the growth of 


The acoustic pulse generated by the initial distur- 
bance propagates at a speed equal to the sound 
speed plus the flow velocity. Thus, the pulse leaves 
the small interior computation domain (50 x50) very 
quickly. Figure 7a shows the pressure contours at 
t = 140. At this time, the acoustic pulse is gone. 
The contours are associated with the excited un- 
stable waves of the acoustic mode. These unsta- 
ble waves move at a slow speed. Figure 7b is at 
t — 200. On comparing figures 7a and 7b, it is ev- 
ident that there is significant growth of the unsta- 
ble waves. Upon reaching the outermost boundary 
of the computation domain the unstable waves are 
reflected back as short waves. This is illustrated 
in figure 7c. The reflected short waves propagate 
at ultrafast speed. They contaminate the computa- 
tion domain in a short period of time as shown in 
figure 7d. Figure 8 shows the growth of the pres- 
sure waveform of the unstable acoustic mode waves 
in the PML before they reach the outer boundary 
of the computation domain. The measured growth 
rate has been found to agree with that calculated by 
the dispersion relation. 

3. Development of a Stable PML 
3.1. Artificial Selective Damping 

To ensure practicality, the thickness of a PML 
would normally be limited to around 15 to 20 mesh 
spacings. For a PML with such a thickness, it is easy 
to show that if the transmitted wave from the com- 
putation domain is to be reduced by a factor of 10 5 
in the presence of a subsonic mean flow, the damp- 
ing coefficient a of (1) should have a value of about 
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1.5. By solving the dispersion relations (2) and (3) 
numerically, it has been found that for a = 1.5 the 
unstable wave solutions have only a modest rate of 
growth. Moreover, these waves, generally, have short 
wavelengths. Mild instabilities of this type can be 
effectively suppressed by the addition of artificial se- 
lective damping terms 12,14 to the discretized govern- 
ing equations. The advantage of using artificial se- 
lective damping is that the damping is confined pri- 
marily to short waves. Thus, the perfectly matched 
condition is not adversely affected for the long waves 
(the physical waves) of the computation. 

Consider the first equation of (1). Let (£, m) be 
the spatial indices in the x- and y-directions. The 
semi-discretized form of this equation using the DRP 
scheme with artificial selective damping terms added 
to the right side is, 

d 3 

3r(ui)*,m + + £ a,j[M x (ui — U2)i+j >m \ 

dt a 


where 

£>(«,/?)= £ djieV +€*”). (15) 

j = - 3 

Contours of the damping function D(ct,/ 3) in the 
a — /?- plane are shown in figure 9. 

To demonstrate that suppression of the unstable 
solutions can be achieved by adding artificial selec- 
tive damping terms to the discretized form of equa- 
tion (1), let us consider the unstable solution with 
growth rate given by figure 3. On combining the 
growth rate of figure 3 and the damping rate of fig- 
ure 9 with Ra = 1.421, the resulting growth con- 
tours are shown in figure 10. Outside the dotted 
lines (wavenumber inside the vertical dotted lines 
corresponds to wavelengths too long to fit into a 15 
mesh spacing PML) the combined effects result in 
damping of the waves. Thus all the instabilities of 
the PML equations are effectively suppressed. 


4* (Pi + P2)l+j ,m] (12) 

1 3 

= --5— £ + (U2)*,i+m] 

R * 3 

where dj y s are the artificial selective damping 
coefficients 14 and R& = Ooo ~ is the artificial mesh 
Reynolds number. Terms similar to those on the 
right side of (12) are to be added to all the other 
discretized equations. 

For the purpose of suppressing unstable solutions 
in the PML, we recommend the use of a damping 
curve with a slightly larger half-width then those 
given in ref. [14]. In this work, the following damp- 
ing coefficients (half- width = 0.357r) are used. 


d 0 = 0.3705630354 


di = d_i = -0.2411788110 


d 2 = d. 2 = 0.0647184823 


(13) 


d 3 = d^z = -0.0088211899. 


The damping rate of the artificial selective damp- 
ing terms can be found by taking the Fourier trans- 
form of the right side of (12) (see [12]). Let (a,/?) 
be the transform variables in the (x,y)-plane. The 
rate of damping for wavenumber (a,/?) is, 


damping rate 


i?A 


D{a y p) 


(14) 


3.2. Distributions of a and R A l in the PML 

In the implementation of PML as an absorbing 
boundary condition, Hu 14 suggested letting a vary 
spatially in the form, 


(16) 

where D is the thickness of the PML, d is the dis- 
tance from the interface with the interior domain 
and A is a constant. With the inclusion of artifi- 
cial selective damping, we have found that the use 
of a well-designed smooth distribution of <7 and R^ 1 
at the interface region is important if the perfectly 
matched condition is to be maintained in the finite 
difference form of the system of equations. 

Figure 11 shows a distribution of a and R^ 1 we 
found to work well with the 7-point stencil DRP 
scheme. The R^ 1 curve is zero for the first two 
mesh points closest to the interface. It attends its 
full value (ii^ 1 )max at the 6 th mesh point. A cubic 
spline curve is used in the transition region. With 
this arrangement, the first point that artificial damp- 
ing occurs is the third point from the interface. This 
allows the use of the 7-point symmetric damping 
stencil in the PML except the last three points at 
the outer boundary. For these points, the 5-point 
and the 3-point stencil 14 should be used instead. 

The <7 curve begins with the value a = 0 at 
the fifth mesh point from the interface. The full 
value ( 7 ma x is reached at 8 mesh points further away. 
Again, a cubic spline curve is used in the transition 
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region. The choice of starting the a curve at the fifth 
point is to ensure that the R^ 1 curve has attained 
its full value when a becomes nonzero. 

3.3. Numerical Examples 


P 


dv _ dv 


+ 


w dv 
r d(j) 


2 ww 


w 2 

~ p ~ 


dp 

dr 


(18b) 


To demonstrate the effectiveness of using artificial 
selective damping terms to suppress the instabilities 
of the PML equations, the numerical examples of 
section 2.3 are reconsidered here. Artificial damping 
is now included in the simulations. Figure 12 shows 
the it-contours of the vorticity waves (M x = 0.3, 
My = 0.2, <r m = 10, (^‘) m ax = 10) as they are 
convected from the interior domain to the PML. The 
vorticity wave packet is steadily damped. No sign 
of unstable waves of the type shown in figure 5 is 
detected. Figure 13 shows the corresponding wave- 
form of u at a few selected times. It is clear that 
the pulse is damped continuously once it propagates 
into the PML. The case of the acoustic disturbance 
has also been repeated with similar results. Based 
on these findings, it is concluded that a stable PML 
can be developed by the inclusion of artificial selec- 
tive damping. Such a PML performs very effectively 
as an absorbing boundary condition in an open do- 
main. 


4. PML in Ducted Environments 

We will now consider the use of PML inside a cir- 
cular duct of radius R . Dimensionless variables with 
respect to length scale i?, velocity scale a t (speed of 
sound at r = /?), time scale density scale p t 
(mean density at r — R) and pressure scale p t a * will 
be used. The velocity components in the (z,r,<£) 
directions of a cylindrical coordinate system are de- 
noted by (u, f , For an inviscid compressible flow, 
the most general mean flow (designated by an over- 
bar) is 

u = u(r), F=0, w = w(r), p = p(r) 


dw _ dw dw w dw w v 
dt + U dx + V dr ^ r d<t>+ r 


1 dp 
r d<j) 


(18c) 


du _ du du w du 

dt U dx + V dr r d<f> 


=-! < i8d > 


dp dp w dp p if 2 

m +u di + 7di + ~ v 


+ 7P 


1 dvr 1 dw du 
r dr + r d<t> + dx J 


= 0 


(18e) 


where 7 is the ratio of specific heats. The boundary 
condition at the duct wall is 


r=l, v = 0. (19) 

Solutions of (18) and (19) representing propagat- 
ing wave modes in the duct may be written in the 
form, 


~ p~ 


/ 

' p( r ) ' 

"j 

u 



u(r) 


V 

= Re < 


v{r) 

exp[i(Arx + m<p — wt)] ) 

w 



w(r) 


-p- 



- P( r ) - 

J 


Substitution of (20) into (18) and (19) leads to the 
following eigenvalue problem. 


~ i d 

P + -T-(P vr ) 

u> r dr 


mw _ k u ~ 

P P 

lot u 


p=- j P^-dr + po. (17) 

r 

Small amplitude disturbances superimposed on 
mean flow (17) are governed by the linearized Eu- 
ler equations. They are, 


_ / m _ k J\ 

— p I — w -1 — u = 0 
\wr u> J 


(21a) 



h _ mw\ ~ . 2 tftf 

— u v —1 

w Ljr J u>r 


ipur 

ujr 


dp 

dt 


Id. x if dp _ dp 
+ - r g-r i7vr) + Tdi + u di 


r d<j) 


i dp 
u> dr 


(21b) 


+ P 


f 1 dw 3u\ 

\r d<j> dx ) 


= 0 


(18a) 


P 


[( 


k _ mw 
1 u 


U) wr 



^ dw i w _ 
V — + — V 

dr wr 
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m _ 
= — P 


(21c) 


_ / k _ mw\ ~ i ~ da k _ 

p 1 u w + = —p (21 d) 

w u;r / warj u; 

( fcu muA _ i pw 2 _ 

1 )p + -- t; 

w wr / w r 


+ 7 P 


i cf(v r) m _ k „ 

: — L w u 

ujt dr ur uf 


r = 1, v = 0 


= 0 (21e) 

( 22 ) 


For a given azimuthal mode number m and fre- 
quency u;, (the wavenumber) is the eigenvalue. 
Corresponding to an eigenvalue is an eigenvector 
[/?, 5, v, w t p \ , which describes the radial profile of the 
wave mode. 

4.1. Perfectly Matched Condition in Ducted 
Flows 

Suppose a perfectly matched layer is to be set 
up as a termination boundary of a computation do- 
main inside a duct. By splitting the variables; e.g., 
p — pi + p2i etc. in the standard manner, the 
PML equations corresponding to the linearized Eu- 
ler equations ((18a) to (18e)) are, 

dpi . 1 d r _, . s , , w d(pi + p 2 ) 


p d(w 1 + W2 ) _ 
r d<£ 


= 0 


(23a) 


dp 2 , , _ d(p i + p 2 ) 

«T + ,, « + “ ai 


w , . 

+ — (vi + t> 2 ) 


1 d(pi + p 2 ) 
r 


90 


9u? 2 _ 9(lt)i + tt»2) 

_ +<tw)2+u __ — 


= 0 


(23e) 

(23f) 


dt*i , , du w _d{u i + tx 2 ) 


= 0 


(23g) 


du 2 , , _ d(m + u 2 ) 

■sr + " J + ” ax'" 


__a(pi_±P2> 

9x 

dpi , wd(p\ +p 2 ) , pw 2 ,.. , .. , 

aT + 7 + ~ ( ' 11 + ” 2) 


(23h) 


+ 7P 


1 3r(vi + v 2 ) 1 9(wi + w 2 ) 

r dr r d<fr 


= 0 (23i) 

dp 2 . . _ d(pi + p 2 ) , _ d{u\ + u 2 ) 

+ 7P 

= 0 (23j) 

where a is the damping coefficient in the PML. The 
boundary condition is 


r = 1, 


v\ 4- V 2 — 0. 


(24) 


In the PML, the duct modes are represented by 
solutions of the form (similar to (20)), 


d(ui + U 2 ) _ „ 

a* -° 


(23b) 

Pi (r, x, t) = Re [p t ( r ) e *(-+ m ^“‘)] , (25) 

\dv\ m w d(v\ + v%) 2 u>(tui + W 2 ) 


etc., where k is the wavenumber. On substituting 
(25) into (23) and (24) and on defining 

^ dt r 90 

r 


w 2 9 , v 

~ (Pi + P 2 ) — = “^(Pi +P 2 ) 

(23c) 

P = Pi + p 2 

U — U\ -\-U2 

du 2 _d(t/i+i> 2 ) 

— + <™ 2 + u 
dt dx 

= 0 

(23d) 

v = v\ +V 2 (26) 

dw\ dw w d(wi + w\) 

-jT + ('i + •>)■* + , at. 


W = Wi + U>2 

P = P1+P2 
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it is straightforward to find that the duct modes in 
the PML are given by the solutions of the following 
eigenvalue problem. 


^ i d . ^ . mw ^ ku 

p+ — -ripvr) p — , 

wr dr u)r w + iu 


( m ^ k _ \ n 

— uH —r- u = 0 

ojr w + iG I 


(27a) 


the interior region of the computation domain is the 
same as that in the PML assures that there is perfect 
matching. That is, a propagating duct mode inci- 
dent on the PML will be totally transmitted into 
the PML without reflection. If the mean flow is 
nonuniform, some of the duct modes may involve 
Kelvin-Helmholtz or other types of flow instability 
waves. However, the perfectly matched condition is 
still valid for these waves. 


( k _ mw \ ^ .2 ww 

1 u v —i 

w + ta wr I wr 


i pw 2 

wr 


i dp 
t j dr 


(27b) 


( k _ mw\ ^ i ^ dw 

1 — u w + — v — 

w 4- icr wr J w dr 


i w ^ 

4 v 

wr 


m ^ 
= — V 


[(■- 


k _ m w \ ^ 

— u u 

w + icr wr I 


i ^du 

+ -V — 
w dr 


W 4“ 1(7 


— p 


(27c) 


(27d) 


f — _\ 2 

ku mtii L i pw ^ 

1 ; P + V 

w 4 - %<j wr I w r 


+ 7 P 


i d(vr) m ^ 

w — 


wr dr wr w + icr 


4.2. The Case of Uniform Mean Flow 


From (25) and (29), the transmitted wave mode 
has the form 


(30) 

If the wave mode is nondispersive, then the in- 
verse of the phase velocity, is positive for waves prop- 
agating in the x-direction and negative in the op- 
posite direction. For these nondispersive waves, the 
transmitted waves are spatially damped; a condition 
needed by the PML if it is to serve as an absorbing 
boundary condition. However, inside a duct, the 
wave modes are dispersive. The direction of prop- 
agation is given by the group velocity We will 
now show that in the presence of a uniform mean 
flow there is a band of acoustic duct modes for which 
the group velocity and the phase velocity have op- 
posite signs. Therefore, for this band of waves, the 
transmitted waves would grow spatially instead of 
being damped. 

By eliminating all the other variables in favor of 
p(r), it is straightforward to find, in the case of a 
uniform mean flow of Mach number M, (21) and (22) 
reduce to the following simple eigenvalue problem. 


d 2 P 1 d£ 

dr 2 r dr 


(w - Mkf - k 


2 



p = 0(31) 


= 0. (27e) 

The boundary condition is 

r = 1, v — 0. (28) 

The eigenvalue is k . On comparing eigenvalue prob- 
lem (21) and (22) with eigenvalue problem (27) and 
(28), it is immediately clear that they are the same if 
£ in (21) is replaced by • Thus the eigenvalues 
are related by 


d p 

(32) 

The eigenfunction is 


P — J m(^mn 

(33) 

where J m ( ) is the m th order Bessel function and 
X mn is the n th root of 

J'Mmn) = 0. 

(34) 


K = Jfe(l + £). (29) 

On the other hand, the eigenvectors are identical. 
The fact that the eigenvectors of a duct mode in 


By substitution of (33) into (31), it is found that 
the dispersion relation or eigenvalue equation for the 
(m,n) th acoustic duct mode is 

(w — Mk) 2 — k 2 — A^ n . (35) 


8 



The axial wavenumber of the mode at frequency cj initial condition is, 
are given by the solution of (35). They are, 


k ± 


-h,M±[h> a -(l-M a )AS, n ]» 

(1 - A/ 2 ) 


(36) 


The group velocity of the duct mode rriay be deter- 
mined by implicit differentiation of (35). This gives, 


d» = ^-(1-M 2 )AJ^(1-M 2 ) 
dk 

In (37), the upper sign corresponds to k = k+ and 
the lower sign corresponds to k — k-> For subsonic 
mean flow, clearly ^ > 0 for k — k+ and < 0 
for Jfc = ifc_. Therefore, the downstream propagat- 
ing waves have wavenumber given by k = while 
the upstream propagating waves have wavenumber 
equal to fc_. 

From (37), it is easy to show that for (1 — 
M 2 )%\ mn < uj < A mn the phase velocity ^ is nega- 
tive although the group velocity is positive. Accord- 
ing to (29), for waves in this frequency band, the 
transmitted wave in the PML will amplify spatially. 
This renders the PML useless as an absorbing layer 
except for M = 0. In the absence of a mean flow 
normal to the PML ( M = 0), Ar + will not be neg- 
ative by (36). Thus, the transmitted waves in the 
PML are evanescent. For this special condition, the 
PML can again be used as an absorbing boundary 
condition. 


4.3. Numerical Examples 

To demonstrate that a PML in a ducted environ- 
ment actually supports a band of amplifying wave 
modes, a series of numerical simulations has been 
carried out. In the simulations, a uniform mesh with 
Ax = A r = 0.04 covering the entire computation 
domain from x = —6.0 to x = 12.0 is used. The 
PML in the upstream direction begins at x = —3.0 
and extends to x = -6.0. In the downstream di- 
rection, the PML occupies the region from x = 3.0 
to x = 12.0. The dimensionless damping constant 
(nondimensionalized by ^ ) cr is set equal to 25.0. 
The results of two simulations, one with a mean flow 
Mach number 0.4, the other with no mean flow are 
reported below. 

For convenience, only the axisymmetric duct 
modes are considered. The computation uses the 
7-point stencil DRP scheme 13 . The acoustic distur- 
bances in the computation domain is initiated by a 
pressure pulse located at x = 0 and r = 0.5. The 


t = 0, u = v = 0, 


p — p ~ exp 


~(*n 2 ) 


(x 2 + (r — 0.5) 2 )1 
16 


( 38 ) 


Figure 14 shows the time evolution of the acous- 
tic disturbance inside the computation domain at 
M = 0.4. Specifically, the pressure waveforms along 
the line r = 0.38 are shown at t = 10, 13, 15 and 16. 
As can be seen, once the pressure pulse is released, 
it spreads out and propagates upstream and down- 
stream. Figure 14a indicates that at time t = 10 the 
front of the acoustic disturbance has just entered the 
PML in the downstream direction. There is no evi- 
dence of wave reflection at the interface between the 
PML and the interior computation domain. The 
transmitted wave grows spatially as shown in fig- 
ure 14b. The amplitude of the transmitted wave in- 
creases steadily as they propagate across the PML. 
This is shown in figures 14c and 14d. When the am- 
plified waves reach the outermost boundary of the 
PML, large amplitude spurious waves are reflected 
back. This quickly contaminates the entire compu- 
tation domain. 

Figure 15 shows the same simulation except that 
there is no mean flow. In the absence of a mean 
flow, the PML acts as an absorbing layer. Figure 15a 
shows the entry of the acoustic pulse into the down- 
stream PML. Figures 15b to 15d show the damping 
of the acoustic pulse in time in the PML. The slowest 
components to decay are the long waves. This is in 
agreement with the analysis of the previous section. 

5. Concluding Remarks 

In this paper, we have shown that the application 
of PML as an absorbing boundary condition for the 
linearized Euler equations works well as long as there 
is no mean flow in the direction normal to the layer. 
For open domain problems, the PML equations, in 
the presence of a subsonic mean flow normal to the 
layer, support unstable solutions. The growth rate of 
the unstable solutions is, however, not large. These 
unstable solutions can, generally, be suppressed by 
the addition of artificial selective damping. In the 
case of a ducted environment, we find that because 
of the highly dispersive nature of the duct modes, a 
band of the transmitted waves in the PML amplifies 
instead of being damped. This seemingly renders the 
PML totally ineffective as an absorbing boundary 
condition. 
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One of the important advantages of using an ab- 
sorbing boundary condition instead of other numer- 
ical boundary treatments is that the boundary of 
the computation domain may be put much closer to 
the source of disturbances. In this way, a smaller 
computation domain may be used in a numerical 
simulation. For open domains, such an absorbing 
boundary condition can be developed by the use of 
PML with artificial selective damping terms. Unfor- 
tunately, the same is not possible for internal ducted 
flow. An effective numerical anechoic termination 
for ducted domains has yet to be developed. 
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Figure 5. Numerical simulation showing the generation and prop- 
agation of unstable vorticity-mode waves in the PML. M x = 0.3, 
M v = 0.2, Ofj\ = 1.0. (a) * = 0, (b) t = 50, (c) t = 90, (d) t = 130. 


Contours of the u velocity component. 0.1, 0.05, 

— 0.01, 0.01, 0.05, 0.1. 



Figure 7. Numerical simulation showing the generation and prop- 
agation of unstable acoustic-mode waves in the PML. 

M x — 0.5, M y = 0.0, a m — 1.5. Contours of pressure. 

(a) t = 140 , p = 10" 4 (b) t = 200 , p = 5. 10' 3 

(c) t = 260 , p = 5. 10“ 2 (d) t = 300 , p = 5. 10" 2 



X 
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Figure 6. Waveforms of u showing the generation of unstable 
vorticity-mode waves excited by vorticity waves convected from 
the interior computation domain to the PML and the subsequent 
contamination of the interior computation domain. M x = 0.3, 
My - 0 . 2 , o m = 1 . 0 . 


Figure 8. Waveforms of pressure along y — 0 showing the gen- 
eration of unstable acoustic-mode waves in the PML excited by 
acoustic disturbances from the computation domain. M x = 0.5, 
My = 0.0, o m — 1.5. 
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Figure 9. Contours of constant D(ar,/?) in the a— (3 plane. Damping 
coefficients dj f s are given by (13). 



-3.0 -2.0 -t.O 0.0 1.0 2.0 3.0 


o 

Figure 10. Contours of combined growth and damping rates. 

M x = 0.3, M v = 0.0, a = 1.5, R& = 1.42. Damping coefficients 
dj’s are given by (13). 



X x 


Figure 12. Damping of a vorticity wave packet in the PML in- 
cluding artificial selective damping terms. M x = 0.3, M v = 0.2, 
a m = 1.0, (Rl')ma: = 1.0. (a) t = 0, (b) t = 50, (c) t = 90, 

(d) t = 170. Contours of the u velocity component. 

0.1, 0.05, 0.01, 0.01, 

0.05, 0.1. 
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Figure 13. Waveforms showing the damping of a vorticity wave Figure 15. Pressure waveforms along the line r = 0.38 of a circular 
packet as it is convected into a PML with artificial selective damp- duct without mean flow showing the damping of an acoustic pulse 
ing terms. M x = 0.3, M v = 0.2, a m = 1.0, = 1.0. in the PML. Ax = Ar = o m = 25.0. 



X 

Figure 14. Pressure waveforms along the line r = 0.38 of a circular 
duct with uniform mean flow at Mach 0.4 showing the excitation 
and growth of the unstable solution in the PML by an acoustic 
pulse. Ax = Ar = o m = 25.0. 


14 






THE EFFECT OF NOZZLE GEOMETRY ON THE 
NOISE OF HIGH-SPEED JETS 


Christopher K.W. Tam 


Department of Mathematics, Florida State University 
Tallahassee, FL 32306-4510 USA 


ABSTRACT 

This paper examines the effectiveness of jet noise reduction by the use of different 
nozzle exit geometry. Since there will be thrust loss associated with a nozzle of complex 
geometry, consideration is confined to practical configurations with reasonably small 
thrust loss. In this study, only jets with a single stream are considered. The nozzle 
configurations examined are circular, elliptic and rectangular. Included also are plug 
nozzles as well as a suppressor nozzle. It is shown that the measured turbulent mixing 
noise of the jets from these nozzles consists of two independent components. The noise 
spectrum of each component is found to fit the shape of a seemingly universal similarity 
spectrum. It is also found that the maximum levels of the fitted noise power spectra of 
the jets are nearly the same. This finding suggests that nozzle geometry modification 
may not be an effective method for jet noise suppression. 

1. INTRODUCTION 

Reducing high-speed jet noise is currently a high priority research and development 
effort of the aircraft industry. Despite many years of jet noise research, noise reduction is 
a highly empirical endeavor. Since the early work of Westly and Lilley 1 many attempts 
have been made to modify the shape of the nozzle exit in the belief that this would 
reduce the turbulence intensity of the jet leading to a reduction in the radiated noise. 
On following this concept, plug nozzles, corrugated nozzles as well as nozzles with multi- 
chute elements have been introduced for noise suppression purpose. 

The objective of this paper is to examine the effectiveness of jet noise reduction 
by nozzle exit geometry modification. Of course, there will be thrust loss in using a 
nozzle with complex geometry. Our consideration is, therefore, confined to practical 
geometries for which the thrust loss is reasonably small. In order to focus attention 



on nozzle geometry alone, we will only consider jets formed by a single stream. Multi- 
stream jets, invariably, would introduce thermodynamic and other flow parameters as 
variables. Under this circumstance, a simple statement on the effectiveness of nozzle 
configuration for noise suppression cannot be easily made. 

In Section 2 of this paper, the effect of nozzle geometry on the turbulent mixing 
processes in jets is discussed. For high-speed jets the mixing process is influenced only 
by upstream events. Thus the normal expectation is that the nozzle exit configuration 
would exert considerable influence on the development of the large and fine scale turbu- 
lence of the jet flow and hence its noise. In Section 3, turbulent mixing noise data from 
a variety of nozzles will be examined and analyzed. It will be shown that the noise level 
is, to a large extent, insensitive to the nozzle shape. This is true even for jets embedded 
in open wind tunnel flows simulating forward flight effects. This result seems to suggest 
that modification of a nozzle exit configuration may not be an effective method for noise 
suppression. 


2- NOZZLE GEOMETRY AS AN INITIAL CONDITION 

Tam and Chen 3 , based on their observation of the noise directivity and spectrum 
measurements of Seiner et a/. 4 , were the first to clearly suggest that turbulent mixing 
noise from high-speed jets is made up of two components. One component is in the form 
of Mach wave radiation generated by the large turbulence structures of the jet flow. 
This component radiates only in the downstream direction. The other component is 
generated by the fine scale turbulence of the jet. The radiated noise has a more uniform 
directivity. Experimental confirmation of the existence of the two noise components 
was not available until the recent investigation of Tam, Golebiowski and Seiner 2 . By 
analyzing the entire data bank of axisymmetric jet noise spectra measured in the Jet 
Noise Laboratory of the NASA Langley Research Center, they were able to extract the 
shapes of two self-similar spectra from the data. They then demonstrated that all the 
noise spectra were made up of a combination of the two similarity spectra. Let S be the 
noise power spectrum ( S has the dimensions of pressure squared per unit frequency) 
then S can be expressed in the following similarity form, 

S ={ AF (-k) +BG (fr)\ (^) 2 (1) 

where F and G 3X6 similarity spectra of the large turbulence structure 

noise and the fine scale turbulence noise respectively, fi is the frequency at the peak 
of the large turbulence structures noise spectrum and /r i s the frequency at the peak 
of the fine scale turbulence noise spectrum. The spectrum functions are normalized 
such that F(l) = G(l) = 1. In equation (1), A and B are the amplitudes of the 
independent spectra. They have the same dimensions as S . Dj is the fully expanded 
jet diameter and r is the distance between the noise measurement point and the nozzle 
exit. The amplitudes A and B and the peak frequencies /l and }f are functions of 
the jet operating parameters and the direction of radiation \ (measured from 

the jet inlet). Vj and a <*> are the jet velocity and the ambient sound speed. T r and Too 
are the reservoir and ambient temperature. One remarkable feature of the similarity 
spectra is that they fit the data well regardless of jet velocity, jet temperature, direction 



of radiation, and whether the jet is perfectly or imperfectly expanded (in the case of 
supersonic jets). These spectra axe used extensively in the present investigation. 

In high-speed jet flows, there is practically very little upstream influence. Thus 
the turbulence level near the end of the core region, where most of the jet noise is 
generated, is affected primarily by the mixing processes upstream and the conditions at 
the nozzle exit. From this point of view, the nozzle geometry may be regarded as an 
initial condition on the spatial evolution of the jet velocity profile and the turbulence 
intensity and spectral content downstream. For noise suppression purposes, the crucial 
question to ask is how sensitive the turbulence level of the jet flow near the end of 
the potential core is to the initial condition at the nozzle exit. There is no question 
that by changing nozzle geometry the entrainment flow and hence jet turbulence in 
the region immediately downstream of the nozzle exit is affected. However, turbulent 
mixing is a highly nonlinear process. It is known, nonlinear process can lead to the 
same asymptotic state regardless of initial conditions. (For a discussion of the lack of 
influence of initial conditions on self-similar turbulent flows, see the work of Tam and 
Chen 5 .) For high Reynolds number jet flows, it is possible that a jet issued from a 
noncircular nozzle evolves quickly into a more or less axisymmetric jet before the end 
of the core is reached. In such a case, the radiated noise would be similar to that of a 
circular jet both in intensity and spectral content. In the next section, it will be shown 
that this appears to be the case. 

3. EVALUATION AND COMPARISONS OF DATA 

Supersonic jet noise data from two sources are used in the present study. The first 
set of data is taken from the data bank of the Jet Noise Laboratory of the NASA Langley 
Research Center. This set of data consists of noise spectra from a Mach 2 aspect ratio 
3 elliptic jet and a Mach 2 aspect ratio 7.6 rectangular jet. These are high quality data; 
comparable to those used in the work of Tam, Golebiowski and Seiner 2 . 

The second set of data is taken from the published measurements of Yamamoto et 
a/. 6 . In this series of experiments, six nozzles are used. They include a conical nozzle, 
a convergent- divergent (C-D) round nozzle, a convergent annular plug nozzle, a C-D 
annular plug nozzle, a 20-chute annular plug suppressor nozzle with convergent flow 
segment terminations and a 20-chute annular plug suppressor nozzle with C-D flow el- 
ement terminations. The noise spectra of the jet from the fifth nozzle, however, are 
strongly different from the same configuration suppressor nozzle but with C-D flow ele- 
ment terminations and the other nozzles. Without knowing the cause of the difference, 
it is decided to ignore the data associated with this nozzle. 

3.1 COMPARISONS WITH SIMILARITY NOISE SPECTRA 

Figure 1 shows direct comparisons between the measured elliptic and rectangular 
jet noise spectra at Mach 2 and = 1.8 from the NASA Langley Research Center 
and the similarity spectrum for the large turbulence structures noise of Tam et a/. 2 at 
X = 150 deg. The elliptic jet noise data are measured on three planes containing the 
jet axis. One is on the minor axis plane, one on a plane at 58 degrees to the minor 
axis plane and the third on the major axis plane. They are the top three curves in the 
figure. The bottom two curves are from the rectangular jet noise data measured on the 
minor and major axis planes. As can be seen, there is good agreement between the 
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measured spectrum shapes and the similarity noise spectrum (the F(-£) function of 
equation (1)). This is so despite the fact that the nozzle geometries are very different. 

Comparisons between the measured spectra at x = 90 deg. and the similarity noise 
spectrum or the fine scale turbulence noise (the G(4^) function of equation (1)) for the 
elliptic and rectangular jets are given in Figure 2. Again, the top three curves are those 
of the elliptic jet and the bottom two curves are of the rectangular jet measured on the 
same azimuthal planes as in Figure 1. It is evident that there is good agreement overall 
regardless of nozzle shapes. 

Figure 3 shows the noise spectrum shapes of the Yamamoto et al. data 6 at x = 150 
deg. The jet velocity in each case is very close to 2420 ft /sec and the total temperature 
is approximately 1715 deg. Rankine. The four spectra are (from the top down) from 
the C-D round nozzle, the convergent annular plug nozzle, the C-D annular plug nozzle 
and the 20-chute annular suppressor nozzle. The data from the conical nozzle is nearly 
the same as the C-D round nozzle and is, therefore, not displayed. The full curves are 
the similarity noise spectrum (the F( -j^) function) of Tam et a/. 2 . On ignoring the very 
low frequency part of the noise spectrum, it is clear that the agreement between the 
measured data and the similarity spectrum is good for all the cases. 

Figure 4 shows similar comparisons as in Figure 3 but at x = 90 deg. By compar- 
ing the several spectra shown, the facility noise contamination at low frequencies can 
be readily detected. The full curves are the similarity spectrum given by the G(-£) 
function. Overall, there is again good fit between the data and the similarity spectrum. 

3.2 COMPARISONS OF MAXIMUM SOUND PRESSURE LEVELS 

To assess whether nozzle geometry has significant influence on high-speed jet noise, 
we compare the sound pressure levels at the peaks of the fitted noise spectra, SPL m&x ^ 
in dB/Hz at r = lOODj from the various jets with the level of the simple circular C-D 
nozzle. The results are shown in Tables 1 to 4. 

Table 1 compares the SPL max of the elliptic jet at temperature ratio (^) of 1.0, 
1.37, 1.80 and 2.27 at jet Mach number 1.98 with the corresponding values of a circular 
jet. We have chosen the microphone measurements at X ” 150 deg. to characterize 
the large turbulence structures noise component and the microphone measurements at 
X = 90 deg. to characterize the fine scale turbulence noise component. The first row 
of data is measured in the minor axis plane. The second row is measured in a plane 
at 58 degrees to the minor axis plane. The third row is measured in the major axis 
plane. The last row is the data from a circular jet at the same jet velocity and total 
temperature. Within experimental uncertainty, it is clear from the table that the noise 
from the elliptic jet is, first of all, quite axisymmetric. Further, it is nearly the same 
as the circular jet. Table 2 provides direct comparisons between the SP Lmax of the 
rectangular jet and a circular jet. Again, within experimental uncertainty, there is very 
little difference in the noise levels. 

Tables 3 and 4 show the SPLmax data at x = 150 and 90 deg. for the various 
nozzles of the Yamamoto et al. experiments. It is worthwhile to remind the readers 
that the data are converted from ^ octave band measurements and possibly slightly 
contaminated by shock and facility noise. The experimental uncertainty could be as 
large as 2 to 3 dB by our estimate. By comparing all the data with those of the C-D 
nozzle, it is evident that the differences are well within the experimental uncertainty. 
Thus, in spite of the large differences in nozzle geometry, the noise from supersonic jets 



are remarkably the same. Based on these results, it is possible to surmise that nozzle 
exit geometry may not have significant control over the noise of high-speed jets. 

4. CONCLUSION 

Extensive comparisons between the noise radiated by supersonic jets operating at 
various temperatures and velocities with and without simulated forward flight and the 
noise from a circular jet at the same conditions have been carried out. Seven nozzles 
of practical geometries are included in the study. It is found that regardless of nozzle 
geometry, turbulent mixing noise of all the jets is comprised of two components. One 
component is the noise from the large turbulence structures and the other is noise 
from the fine scale turbulence of the jet flow. Further, the radiated sound is largely 
axisymmetric and that the shapes of the spectra of the two noise components are nearly 
the same as those of the similarity spectra of Tam, Golebiowski and Seiner 2 . In addition, 
the noise levels axe essentially independent of nozzle configuration. Based on these 
results, it is concluded (bearing in mind the limited scope of this study) that nozzle 
geometry modification may not be an effective method for jet noise suppression. 
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Table 1. Elliptic jet (aspect ratio 3, Mj = 1.98) 



v = 90 dee. 

y = 141 de 
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Tr/Too 

1.00 

1.37 

1.80 

2.27 

1.00 

1.37 

1.80 

2.27 

measurement plane 

S P L ma .x 

at r = 100-Dj 
(dB/Hz) 

74.3 

75.5 

77.0 


96.8 

99.5 

101.7 


minor axis plane 

74.3 

75.7 

76.8 

78.3 

96.1 

98.8 

100.3 

101.3 

58 deg. plane 

74.5 

75.5 

77.0 

78.6 

94.4 

97.5 

101.7 

101.7 

major axis plane 

75.5 

76.2 

77.3 

78.5 

97.3 

99.3 

100.7 

102.1 

circular jet 


Table 2. Rectangular jet (aspect ratio 7.6, Mj = 2.0) 
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1.82 

2.26 

1.10 

1.82 
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measurement plane 
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at r = 100-Dj 
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76.9 
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minor axis plane 

74.9 

75.9 

77.0 

98.1 

100.2 

100.6 

major axis plane 

76.0 

77.7 

78.8 

98.4 

101.5 

102.6 

circular jet 
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Table 3. Yamamoto et al. data 
( vj ~ 2420 ft/sec, T r ~ 1715 deg. R) 


nozzle 

conical 

C-D nozzle 

convergent 

C-D plug 

suppressor 

inlet angle 

type 

nozzle 

M d = 1.4 

plug nozzle 

nozzle 

nozzle 

X, degree 

SP L -[ nax a t 

98.8 

97.7 

98.7 

99.0 

97.4 

150 

r = 100Dj (dB/Hz) 

77.6 

75.0 

76.6 

77.2 

74.5 

90 


Table 4. Yamamoto et al. data 
( vj ~ 1720 ft/sec, T r ~ 870 deg. R) 


nozzle 

C-D nozzle 

convergent 

C-D plug 

suppressor 

inlet angle 

type 

M d = 1.4 

plug nozzle 

nozzle 

nozzle 

X, degree 

SP imax a t 

95.0 

96.2 ^ 

97.1 

92.5 

150 

r = lOODj (dB/Hz) 

70.3 

73.0 

74.0 

70.0 

90 
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Figure 1. Comparisons between elliptic and rectangular jet noise 
data and the similarity spectrum at \ = 150 deg., ^ = 1.8 

Aspect ratio 3 elliptic jet: (a) minor axis plane, (b) 58 degree 
plane, (c) major axis plane. 

Aspect ratio 7.6 rectangular jet: (d) minor, (e) major axis plane. 
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Figure 2. Comparisons between elliptic and rectangular jet noise 
data and the similarity spectrum at \ = 90deg., ^ = 1.8 

Aspect ratio 3 elliptic jet: (a) minor axis plane, (b) 58 degree 
pla ne, (c) major axis plane. 

Aspect ratio 7.6 rectangular jet: (d) minor, (e) major axis plane. 
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Figure 3. Comparisons between Yamamoto et al. data and the sim- 
ilarity spectrum. V} ~ 2420 ft/sec, T r ~ 1715 deg R, \ = 150 deg; 
o data, similarity spectrum, (a) C-D nozzle, (b) conver- 

gent plug nozzle, (c) C-D plug nozzle, (d) 20-chute C-D suppressor 
nozzle. 
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Figure 4. Comparisons between Yamamoto et al. data and the 
similarity spectrum. Vj ~ 2420 ft/sec, T r — 1715 deg R, x = 90 

deg; o data, similarity spectrum, (a) C-D nozzle, (b) conve 

rgent plug nozzle, (c) C-D plug nozzle, (d) 20-chute C-D suppressor 
nozzle. 





